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The electronic ground state of a periodic system is usually described in terms of ex- 
^v^j tended Bloch orbitals, but an alternative representation in terms of localized "Wannier 

functions" was introduced by Gregory Wannier in 1937. The connection between the 
I Bloch and Wannier representations is realized by families of transformations in a con- 

tinuous space of unitary matrices, carrying a large degree of arbitrariness. Since 1997, 
methods have been developed that allow one to iteratively transform the extended Bloch 
orbitals of a first-principles calculation into a unique set of maximally localized Wannier 
functions, accomplishing the solid-state equivalent of constructing localized molecular 
^s^J orbitals, or "Boys orbitals" as previously known from the chemistry literature. These 

developments are reviewed here, and a survey of the applications of these methods is pre- 
sented. This latter includes a description of their use in analyzing the nature of chemical 
bonding, or as a local probe of phenomena related to electric polarization and orbital 
. . magnetization. Wannier interpolation schemes are also reviewed, by which quantities 

^ computed on a coarse reciprocal-space mesh can be used to interpolate onto much finer 

meshes at low cost, and applications in which Wannier functions are used as efficient 
basis functions are discussed. Finally the construction and use of Wannier functions 
i outside the context of electronic-structure theory is presented, for cases that include 

phonon excitations, photonic crystals, and cold-atom optical lattices. 
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I. INTRODUCTION 

In the independent-particle approximation, the elec- 
tronic ground state of a system is determined by spec- 
ifying a set of one-particle orbitals and their occupa- 
tions. For the case of periodic system, these one-particle 
orbitals are normally taken to be the Bloch functions 
VVik(r) that are labeled, according to Bloch's theorem, by 
a crystal momentum k lying inside the Brillouin zone and 
a band index n. Although this choice is by far the most 
widely used in electronic-structure calculations, alterna- 
tive representations are possible. In particular, to ar- 



rive at the Wannier representation (des Cloizeaux 1963 



E3 Kohn, 1959 Wannier 1937), one carries out a unitary 



transformation from the Bloch functions to a set of local- 
ized "Wannier functions" (WFs) labeled by a cell index 
R and a band-like index n, such that in a crystal the 
WFs at different R are translational images of one an- 
other. Unlike Bloch functions, WFs are not eigenstates 
of the Hamiltonian; in selecting them, one trades off lo- 
calization in energy for localization in space. 

In the earlier solid-state theory literature, WFs were 
typically introduced in order to carry out some formal 
derivation - for example, of the effective-mass treatment 
of electron dynamics, or of an effective spin Hamiltonian 
- but actual calculations of the WFs were rarely per- 
formed. The history is rather different in the chemistry 
literature, where "localized molecular orbitals" (LMO's) 



(Boys 


1960 


1966 


Edmiston and Ruedenberg| 1963| Fos- 


ter and Boys 


1960a b 


1 have played a significant role in 



computational chemistry since its early days. Chemists 
have emphasized that such a representation can provide 
an insightful picture of the nature of the chemical bond 
in a material — otherwise missing from the picture of 
extended eigenstates — or can serve as a compact basis 
set for high-accuracy calculations. 

The actual implementation of Wannier's vision in 
the context of first-principles electronic-structure calcu- 
lations, such as those carried out in the Kohn-Sham 



framework of density- functional theory ( Kohn and Sham 



1965), has instead been slower to unfold. A major rea- 
son for this is that WFs are strongly non-unique. This 
is a consequence of the phase indeterminacy that Bloch 
orbitals ?/> n k have at every wavevector k - or, more gener- 
ally, the "gauge" indeterminacy associated with the free- 
dom to apply any arbitrary unitary transformation to the 
occupied Bloch states at each k. This second indetermi- 
nacy is all the more troublesome in the common case 
of degeneracy for the occupied bands at certain high- 
symmetry points in the Brillouin zone, making a par- 



3 



tition into separate "bands", that could separately be 
transformed in Wannier functions, problematic. There- 
fore, even before one could attempt to compute the WFs 
for a given material, one had first to resolve the question 
of which states to use to compute WFs. 

An important development in this regard was the intro- 



duction by Marzari and Vanderbilt (1997) of a "maximal 



localization" criterion for identifying a unique set of WFs 
for a given crystalline insulator. The approach is simi- 
lar in spirit to the construction of LMO's in chemistry, 
but its implementation in the solid-state context required 
significant developments, due to the ill-conditioned na- 



ture of the position operator in periodic systems (Nen- 



ciu 



1991 ), that was clarified in the context of the "mod- 



ern theory" of polarization ( King-Smith and Vanderbilt 



1993 Resta 1994 1. Marzari and Vanderbilt showed that 



the minimization of a localization functional correspond- 
ing to the sum of the second- moment spread of each Wan- 
nier charge density about its own center of charge was 
both formally attractive and computationally tractable. 



In a related development, Souza et al. (2001) generalized 



the method to handle the case in which one wants to 
construct a set of WFs that spans a subspace containing, 
e.g., the partially occupied bands of a metal. 

These developments touched off a transformational 
shift in which the computational electronic-structure 
community started constructing maximally-localized 
WFs (MLWFs) explicitly and using these for different 
purposes. The reasons are manifold: First, WFs, akin 
to LMO's in molecules, provide an insightful chemical 
analysis of the nature of bonding, and its evolution dur- 
ing, say, a chemical reaction. As such, they have become 
an established tool in the post-processing of electronic- 
structure calculations. More interestingly, there are for- 
mal connections between the centers of charge of the 
WFs and the Berry phases of the Bloch functions as 
they are carried around the Brillouin zone. This con- 
nection is embodied in the microscopic modern theory 
of polarization, alluded to above, and has led to impor- 
tant advances in the characterization and understand- 
ing of dielectric response and polarization in materials. 
Of broader interest to the entire condensed matter com- 
munity is the use of WFs in the construction of model 
Hamiltonians for, e.g., correlated-electron and magnetic 
systems. An alternative use of WFs as localized, transfer- 
able building blocks has taken place in the theory of bal- 
listic (Landauer) transport, where Green's functions and 
self-energies can be constructed effectively in a Wannier 
basis, or that of first-principles tight-binding Hamilto- 
nians, where chemically-accurate Hamiltonians are con- 
structed directly on the Wannier basis, rather than fitted 
or inferred from macroscopic considerations. Finally, the 
ideas that were developed for electronic WFs have also 
seen application in very different contexts. For example, 
MLWF's have been used in the theoretical analysis of 
phonons, photonic crystals, cold atom lattices, and the 



local dielectric responses of insulators. 

Here we review these developments. We first intro- 
duce the transformation from Bloch functions to WFs in 
Sec. [TTJ discussing their gauge freedom and the methods 
developed for constructing WFs through projection or 
maximal localization. A "disentangling procedure" for 
constructing WFs for a non- isolated set of bands (e.g., 
in metals) is also described. In Sec. Ill we discuss vari- 
ants of these procedures in which different localization 
criteria or different algorithms are used, and discuss the 
relationship to "downfolding" and linear-scaling meth- 
ods. Sec. |IV| describes how the calculation of WFs has 
proved to be a useful tool for analyzing the nature of the 
chemical bonding in crystalline, amorphous, and defec- 
tive systems. Of particular importance is the ability to 
use WFs as a local probe of electric polarization, as de- 
scribed in Sec.fV] There we also discuss how the Wannier 
representation has been useful in describing orbital mag- 
netization, NMR chemical shifts, orbital magnetoclcctric 
responses, and topological insulators. Sec. |VI| describes 
Wannier interpolation schemes, by which quantities com- 
puted on a relatively coarse k-space mesh can be used 
to interpolate faithfully onto an arbitrarily fine k-space 



mesh at relatively low cost. In Sec. VII we discuss appli- 
cations in which the WFs are used as an efficient basis 
for the calculations of quantum transport properties, the 
derivation of semiempirical potentials, and for describing 
strongly-correlated systems. Sec. VIII contains a brief 



discussion of the construction and use of WFs in con- 
texts other than electronic-structure theory, including for 
phonons in ordinary crystals, photonic crystals, and cold 
atoms in optical lattices. Finally, Sec. |IX| provides a short 
summary and conclusions. 



II. REVIEW OF BASIC THEORY 

A. Bloch functions and Wannier functions 

Electronic structure calculations are often carried out 
using periodic boundary conditions. This is the most 
natural choice for the study of perfect crystals, and also 
applies to the common use of periodic supercells for the 
study of non-periodic systems such as liquids, interfaces, 
and defects. The one-particle effective Hamiltonian H 
then commutes with the lattice-translation operator Tr, 
allowing one to choose as common eigenstates the Bloch 
orbitals | VVik ) : 



[H, Tr] = 



V'nk(r) = u n k(r)e 



ik-r 



(i) 



where M n k(r) has the periodicity of the Hamiltonian. 

Several Bloch functions are sketched on the left-hand 
side of Fig. [T] for a toy model in which the band of inter- 
est is composed of p-like orbitals centered on each atom. 
For the moment, we suppose that this band is an iso- 
lated band, i.e., it remains separated by a gap from the 
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Bloch functions 
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Wannier functions 

w (x) 




w/x) 



FIG. 1 (Color online) Left: Bloch functions associated with 
a single band in ID, at three different values of wavevector k. 
Right: WFs associated with the same band, forming periodic 
images of one another. Blue dots indicate atoms; green curves 
indicate envelopes e lkx of the Bloch functions. Bloch and WFs 
span the same Hilbert space. 



bands below and above at all k. Since Bloch functions 
at different k have different envelope functions e r , one 
can expect to be able to build a localized "wave packet" 
by superposing Bloch functions of different k. To get a 
very localized wave packet in real space, we need to use 
a very broad superposition in k space. But k lives in the 
periodic Brillouin zone, so the best we can do is to choose 
equal amplitudes all across the Brillouin zone. Thus, we 
can construct 



w Q (r) 



V 



(27T) 



dk ^nk(r) 



(2) 



BZ 



where V is the real-space primitive cell volume and the 
integral is carried over the Brillouin zone (BZ). (See 
Sec. II. A. 3 for normalization conventions.) Equation ^ 
can be interpreted as the WF located in the home unit 
cell, as sketched in the top-right panel of Fig. [T] 



More generally, we can insert a phase factor e~ lk R into 
the integrand of Eq. Q, where R is a real-space lattice 
vector; this has the effect of translating the real-space 
WF by R, generating additional WFs such as Wi and 
W2 sketched in Fig. [T] Switching to the Dirac bra-ket 
notation and introducing the notation that Rn refers to 
the WF w n R in cell R associated with band n, WFs can 



be constructed according to ( Wannier 1937 ) 
V 



Rn 



dk e 



-ik-R 



(3) 



(2tt)3 Jbz 

It is easily shown that the | Rn ) form an orthonormal 
set (see Sec. |II.A.3[) and that two WFs | Rn ) and | R'n ) 



vector R R' ( |Blount[ [l962| . Eq. ^ takes the form of 
a Fourier transform, and its inverse transform is 



E 

R 



,ik-R 



Rn) 



(4) 



(see Sec. II. A. 3). Any of the Bloch functions on the left 
side of Fig. [Tj can thus be built up by linearly superposing 
the WFs shown on the right side, when the appropriate 
phases e lkR are used. 

The transformations of Eqs. ([3| and Q constitute 
a unitary transformation between Bloch and Wannier 
states. Thus, both sets of states provide an equally valid 
description of the band subspace, even if the WFs are 
not Hamiltonian eigenstates. For example, the charge 
density obtained by summing the squares of the Bloch 
functions \ip n k} or the WFs |Rn) is identical; a similar 
reasoning applies to the trace of any one-particle opera- 
tor. The equivalence between the Bloch and the Wannier 
representations can also be made manifest by expressing 
the band projection operator P in both representations, 
i.e., as 



P 



V 



(2tt) s 



dk \ij>nk)(il>r> 



BZ 



E 

R 



Rn)(Rn| 



(5) 



WFs thus provide an attractive option for represent- 
ing the space spanned by a Bloch band in a crystal, being 
localized while still carrying the same information con- 
tained in the Bloch functions. 



1. Gauge freedom 

However, the theory of WFs is made more complex 
by the presence of a "gauge freedom" that exists in the 
definition of the ip n k- In fact, we can replace 



or equivalently, 



I ^nk) 



_ giv^Ck) 



4>nk) 



(6) 



(7) 



transform into each other under translation by the lattice 



without changing the physical description of the system, 
with ip n (k) being any real function that is periodic in 
reciprocal space j]] A smooth gauge could, e.g., be denned 
such that Vk|w.nk) is well defined at all k. Henceforth 
we shall assume that the Bloch functions on the right- 
hand side of Eq. ([3| belong to a smooth gauge, since 



1 More precisely, the condition is that ip n (k+G) = i/? n (k) + G-AR 
for any reciprocal-lattice translation G, where AR is a real-space 
lattice vector. This allows for the possibility that ip n may shift 
by 2-7T times an integer upon translation by G; the vector AR 
expresses the corresponding shift in the position of the resulting 
WF. 
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we would not get well-localized WFs on the left-hand 
side otherwise. This is typical of Fourier transforms: the 
smoother the reciprocal-space object, the more localized 
the resulting real-space object, and vice versa. 

One way to see this explicitly is to consider the R = 
home-cell w n o(r) evaluated at a distant point r; using 
Eq. in Eq. m, this is given by J BZ u nk (r)e* k r dk, 
which will be small due to cancellations arising from the 
rapid variation of the exponential factor, provided that 
u n k is a smooth function of k ( Blount 



1962) 



It is important to realize that the gauge freedom of 
Eqs. ^ and ^ propagates into the WFs. That is, 
different choices of smooth gauge correspond to differ- 
ent sets of WFs having in general different shapes and 
spreads. In this sense, the WFs are "more non- unique" 
than the Bloch functions, which only acquire a phase 
change. We also emphasize that there is no "preferred 
gauge" assigned by the Schrodinger equation to the Bloch 
orbitals. Thus, the non-uniqueness of the WFs resulting 
from Eq. Q is unavoidable. 

2. Multiband case 

Before discussing how this non-uniqueness might be 
resolved, we first relax the condition that band n be a 
single isolated band, and consider instead a manifold of 
J bands that remain separated with respect to any lower 
or higher bands outside the manifold. Internal degen- 
eracies and crossings among the J bands may occur in 
general. In the simplest case this manifold corresponds 
to the occupied bands of an insulator, but more generally 
it consists of any set of bands that is separated by a gap 
from both lower and higher bands everywhere in the Bril- 
louin zone. Traces over this band manifold are invariant 
with respect to any unitary transformation among the J 
occupied Bloch orbitals at a given wavevector, so it is 
natural to generalize the notion of a "gauge transforma- 
tion" to 



E 



c4 k 2 1 i>v 



(8) 



Here Uml is a unitary matrix of dimension J that is 
periodic in k, with Eq. Q corresponding to the special 
case of a diagonal U matrix. It follows that the projection 
operator onto this band manifold at wavevector k, 



-Pk = IV'nkKV'nkl = |'0nk)(V'nk| 



(9) 



n=l 



n=l 



is invariant, even though the |-0«k) resulting from Eq 
are no longer generally eigenstates of H, and n is no 
longer a band index in the usual sense. 

Our goal is again to construct WFs out of these trans- 
formed Bloch functions using Eq. (|]). Figs.^a) and (b) 



FIG. 2 (Color online) MLWFs constructed from the four va- 
lence bands of Si (left) and GaAs (right; Ga at upper right, 
As at lower left), having the character of <r-bonded combina- 
tions of sp 3 hybrids. The red and blue isosurfaces correspond 
to two opposite values for the amplitudes of the real- valued 
MLWFs. 



show, for example, what the result might eventually look 
like for the case of the four occupied valence bands of Si 
or GaAs, respectively. From these four bands, one ob- 
tains four equivalent WFs per unit cell, each localized on 
one of the four nearest-neighbor Si-Si or Ga-As bonds. 
The presence of a bond-centered inversion symmetry for 
Si, but not GaAs, is clearly reflected in the shapes of the 
WFs. 

Once again, we emphasize that the gauge freedom ex- 
pressed in Eq. ^ implies that the WFs are strongly non- 
unique. This is illustrated in Fig. [3j which shows an al- 
ternative construction of WFs for GaAs. The WF on the 
left was constructed from the lowest valence band n=l, 
while the one on the right is one of three constructed from 
bands n=2-4. The former has primarily As s character 
and the latter has primarily As p character, although 
both (and especially the latter) contain some Ga s and 
p character as well. The WFs of Figs. [5Jb) and Fig. [3] 
are related to each other by a certain manifold of 4x4 
unitary matrices U^m relating their Bloch transforms in 
the manner of Eq. Q . 

However, before we can arrive at well-localized WFs 
like those shown in Figs.[2]and[3j we again have to address 
questions of smoothness of the gauge choice expressed in 
Eq. ((8). This issue is even more profound in the present 
multiband case, since this smoothness criterion is gen- 
erally incompatible with the usual construction of Bloch 
functions. That is, if we simply insert the usual Bloch 
functions \ipnk), defined to be eigenstates of H, into the 
right-hand side of Eq. ^3]), it is generally not possible to 




X K 



FIG. 3 (Color online) MLWFs constructed from the s band 
(left) or from the three p bands (right) of GaAs. 
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produce well-localized WFs. The problem arises when 
there are degeneracies among the bands in question at 
certain locations in the Brillouin zone. Consider, for ex- 
ample, what happens if we try to construct a single WF 
from the highest occupied band n — 4 in GaAs. This 
would be doomed to failure, since this band becomes de- 
generate with bands two and three at the zone center, 
r, as shown in Fig. || As a result, band four is non- 
analytic in k in the vicinity of T. The Fourier transform 
of Eq. ([3]) would then result in a poorly localized object 
having power-law tails in real space. 

In such cases, therefore, the extra unitary mixing ex- 
pressed in Eq. ^ is mandatory, even if it may be op- 
tional in the case of a set of discrete bands that do not 
touch anywhere in the BZ. So, generally speaking, our 
procedure must be that we start from a set of Hamilto- 
nian eigenstates \tp n u) that are not per se smooth in k, 
and introduce unitary rotations Umn that "cancel out" 
the discontinuities in such a way that smoothness is re- 
stored, i.e., that the resulting |"0 n k) of Eq. ^ obey the 
smoothness condition that VklV'nk) remains regular at 
all k. Then, when these are inserted into Eq. ^ 

in place of the \ip n k), well-localized WFs should result. 
Explicitly, this results in WFs constructed according to 



\Rn) 



V 



(27T) 



dke 



-ik-R 



BZ 



m—1 



t4 k 2i^ mk ) . (10) 



The question remains how to choose the unitary rota- 
tions Umn so as to accomplish this task. We will see that 
one way to do this is to use a projection technique, as 
outlined in the next section. Ideally, however, we would 
like the construction to result in WFs that are "maxi- 
mally localized" according to some criterion. Methods 
for accomplishing this are discussed in Sec. |II.C| 



be the number of unit cells in the periodic supercell, or 
equivalently, the number of mesh points in the BZ, it is 
possible to keep the conventions close to the continuous 
case by defining the Fourier transform pair as 



|V«k> = Yl eik ' R l R "> 



R 



R 



n) = -Y 



(12) 



-ik-R 



\4>nk) 



with (V'nklV'mk') = N5 nm <5kk' j so that Eq. ^ becomes, 
after generalizing to the multiband case, 



P = Jf E 1^*) Wnkl = E WW 
nk nR 

Another commonly used convention is to write 



(13) 



IV>nk) = -j= 

V R 



e ik R |Rra) 



(14) 



R, 



i > = 7^ E e ~ ik R 1 w 



with (i)nk\ipmk') = Snm <W and Eq. (|13) replaced by 

/ 3 = ^|^nk)(^ lk | = ^|Rn)(Rn| . (15) 



nk 



nR 



In either case, it is convenient to keep the |u„k) functions 
normalized to the unit cell, with inner products involving 
them, such as (u TO k|it„k) , understood as integrals over 



one unit cell. In the case of Eq. (14), this means that 



3. Normalization conventions 

In the above equations, formulated for continuous k, 
we have adopted the convention that Bloch functions are 
normalized to one unit cell, f v dr |Vvik(r)| 2 = 1, even 
though they extend over the entire crystal. We also de- 
fine (f\g) as the integral of f*g over all space. With 
this notation, (VViklVVik) is not unity; instead, it diverges 
according to the rule 

(27T) 3 

V 

With these conventions it is easy to check that the WFs 
in Eqs. (3|4) are properly normalized, i.e., (RnjR'm) = 



mk' 



S nm £ 3 (k-k') 



(ii) 



<5rr' S n ... 

It is often more convenient to work on a discrete uni- 
form k mesh instead of continuous k spacej^] Letting N 



B. Wannier functions via projection 

A simple yet often effective approach for constructing 
a smooth gauge in k, and a corresponding set of well- 
localized WFs, is by projection - an approach that finds 



its roots in the analysis of des Cloizeaux (1964a). Here, 



as discussed, e.g., in Sec. IV.G.l of Marzari and Van- 



derbilt ( 1997 ) , one starts from a set of J localized trial 



The discretization of k-space amounts to imposing periodic 



orbitals g n { v ) corresponding to some rough guess for the 



boundary conditions on the Bloch wavefunctions over a super- 
cell in real spa ce. Thus, it s hould be kept in mind that the WFs 
given by Eqs. ( ] 1 2[ i and | |14| l are not truly localized, as they also 
display the supercell periodicity (and are normalized to a super- 
cell volume) . Under these circumstances the notion of "Wannier 
localization" refers to localization within one supercell, which is 
meaningful for supercells chosen large enough to ensure negligible 
overlap between a WF and its periodic images. 
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WFs in the home unit cell. Returning to the continuous- 
k formulation, these <? n (r) are projected onto the Bloch 
manifold at wavevector k to obtain 



\4>nV 



J 

E 



\4>mk)(lpmk\g n ) 



(16) 



which are typically smooth in k-space, albeit not or- 
thonormal. (The integral in (VVnklffn) is over all space as 
usual.) We note that in actual practice such projection 
is achieved by first computing a matrix of inner products 
(Akjmn = (ip m k\g n ) and then using these in Eq. (|l~6|. 
The overlap matrix (5 k ) mn = (0 m k|<?Wv = { A k A k)mn 
(where subscript V denotes an integral over one cell) 
is then computed and used to construct the Lowdin- 
orthonormalized Bloch-like states 



\lpn 



J 

E 



10, 



(17) 



These \ip n k) have now a smooth gauge in k, are related 
to the original \ipnk) by a unitary transformation^] and 
when substituted into Eq. ([3| in place of the \ip n k) re- 
sult in a set of well-localized WFs. We note in passing 
that the jVVsk) are uniquely defined by the trial orbitals 
g n (r) and the chosen (isolated) manifold, since any arbi- 
trary unitary rotation among the IVVik) orbitals cancels 
out exactly and does not affect the outcome of Eq. [16J 
eliminating thus any gauge freedom. 

We emphasize that the trial functions do not need to 
resemble the target WFs closely; it is often sufficient to 
choose simple analytic functions (e.g., spherical harmon- 
ics times Gaussians) provided they are roughly located 
on sites where WFs are expected to be centered and have 
appropriate angular character. The method is successful 
as long as the inner-product matrix does not become 
singular (or nearly so) for any k, which can be ensured by 
checking that the ratio of maximum and minimum val- 
ues of det Sk in the Brillouin zone does not become too 
large. For example, spherical (s-like) Gaussians located 
at the bond centers will suffice for the construction of 
well-localized WFs, akin to those shown in Fig. [2j span- 
ning the four occupied valence bands of semiconductors 
such as Si and GaAs. 



C. Maximally localized Wannier functions 

The projection method described in the previous sub- 



section has been used by many authors ( Ku et al. 2002 



3 One can prove that this transformation is unitary by performing 
the singular value decomposition A = ZDW\ with Z and W 
unitary and D real and diagonal. It follows that A(A^ A) -1 / 2 is 
equal to ZW^ , and thus unitary. 



Lu et al. 


2004 


Qian et al. 


2008 



Stephan et al. 20001, 



as has a related method involving downfolding of the 



band structure onto a minimal basis ( Andersen and Saha- 



Dasgupta 2000 Zurek et al. 2005); some of these ap- 



proaches will also be discussed in Sec. |III.B| Other au- 
thors have made use of symmetry considerations, analyt- 



icity requirements, and variational procedures (Smirnov 



and Usvyat 2001 Sporkmann and Bross 1994 1997). 



A very general and now widely used approach, how- 
ever, has been that developed by Ma rzari and Vanderbilt| 
( 1997 ) in which localization is enforced by introducing a 



well-defined localization criterion, and then refining the 

fk) 

Umn in order to satisfy that criterion. We first give an 
overview of this approach, and then provide details in the 
following subsections. 

First, the localization functional 

n = [( On | r 2 | On ) - ( On | r | On ) 2 ] 

n 



is defined, measuring of the sum of the quadratic spreads 
of the J WFs in the home unit cell around their cen- 
ters. This turns out to be the solid-state equivalent of the 



1960 



Foster-Boys criterion of quantum chemistry (Boys 
1966| |Foster and BoysJ | 1960a|b ). The next step is to 



express Q in terms of the Bloch functions. This requires 
some care, since expectation values of the position opera- 
tor are not well defined in the Bloch representation. The 



needed formalism will be discussed briefly in Sec. II. Gl 
and more extensively in Sec. |V.A~T] with much of the con- 
ceptual work stemming from the earlier development the 



modern theory of polarization ( Blount[ |1962[ |King-Smith 



and Vanderbilt 


1993 Resta 1992, 


1994; Vanderbilt and 


King-Smith 


1993 


1. 



Once a k-space expression for il has been derived, max- 
imally localized WFs are obtained by min imizing it with 

(k) J I 

respect to the U mn appearing in Eq. (10). This is done 



as a post-processing step after a conventional electronic- 
structure calculation has been self-consistently converged 
and a set of ground-state Bloch orbitals | ipmk ) has been 
chosen once and for all. The Umn are then iteratively re- 
fined in a direct minimization procedure of the localiza- 
tion functional that is outlined in Sec. III. Dl below. This 
procedure also provides the expectation values (r 2 ) n and 
f n ; the latter, in particular, are the primary quantities 
needed to compute many of the properties, such as the 
electronic polarization, discussed in Scc.lVj If desired, the 

(k) — 

resulting Umn can also be used to construct explicitly the 
maximally localized WFs via Eq. (10). This step is typ- 
ically only necessary, however, if one wants to visualize 
the resulting WFs or to use them as basis functions for 
some subsequent analysis. 
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1. Real-space representation 



An interesting consequence stemming from the choice 



of ( 18 ) as the localization functional is that it allows 



a natural decomposition of the functional into gauge- 
invariant and gauge-dependent parts. That is, we can 
write 



where 



; On I r 2 | On ) - I (Rm|r|On) | 



R,?n 



and 



n = E |( R H r l° r 

ri Rm/On 



(19) 



(20) 



(21) 



It can be shown that not only ft but also Dt is pos- 
itive definite, and moreover that f^i is gauge-invariant, 
i.e., invariant under any arbitrary unitary transformation 



( 10 1 of the Bloch orbitals ( Marzari and Vanderbilt 1997 ) . 



This follows straightforwardly from recasting Eq. ( 20 1 in 



terms of the band-group projection operator P, as de- 



fined in Eq. (15), and its complement Q — 1 — P: 



fii = ^{Qn\r a Qr a \0i 

not 



(22) 



The subscript 'c' indicates trace per unit cell. Clearly 
f2i is gauge invariant, since it is expressed in terms of 
projection operators that arc unaffected by any gauge 
transformation. It can also be seen to be positive definite 
by using the idempotency of P and Q to write fli = 
E«Tr c [{Pr a Q)(Pr a Qn = T, a \\Pr a Q\\l 

The minimization procedure of Q thus actually cor- 
responds to the minimization of the non-invariant part 
fl only. At the minimum, the off-diagonal elements 
| (Rm|r|0n) | 2 are as small as possible, realizing the best 
compromise in the simultaneous diagonalization, within 
the subspace of the Bloch bands considered, of the three 
position operators x, y and z, which do not in general 
commute when projected onto this space. 



2. Reciprocal-space representation 



As shown by Blount (19621, matrix elements of the 



position operator between WFs take the form 

(Rn\r\0m) = i J dk e^ R (u nk \\/ k \u mk ) (23) 



and 



(Rn\r 2 \0m) = -7^3 J dk e 4k R (u nk |V£|u mk ) . (24) 



These expressions provide the needed connection with 
our underlying Bloch formalism, since they allow to ex- 
press the localization functional £1 in terms of the ma- 
trix elements of V k and V 2 ,. In addition, they allow to 
calculate the effects on the localization of any unitary 
transformation of the |it nk ) without having to recalcu- 
late expensive (especially when plane- wave basis sets are 
used) scalar products. We thus determine the Bloch or- 
bitals |u mk ) on a regular mesh of k-points, and use finite 
differences to evaluate the above derivatives. In particu- 
lar, we make the assumption that the BZ has been dis- 
cretized into a uniform Monkhorst-Pack mesh, and the 
Bloch orbitals have been determined on that mesh0 

For any /(k) that is a smooth function of k, it can 
be shown that its gradient can be expressed by finite 
differences as 

V/(k) = ™b b [/(k + b) - /(k)] + 0{b 2 ) (25) 

b 

calculated on stars ("shells") of near-neighbor k-points; 
here b is a vector connecting a k-point to one of its neigh- 
bors, Wb is an appropriate geometric factor that depends 
on the number of points in the star and its geometry 



(see Appendix B in Marzari and Vanderbilt (1997) and 
Mostofi et al. (20081 for a detailed description). In a 
similar way, 



|V/(k)| : 



w b [/(k + b)-/(k)] 2 +0(6 3 



(26) 



It now becomes straightforward to calculate the ma- 
trix elements in Eqs. (23) and (24). All the information 



needed for the reciprocal-space derivatives is encoded in 
the overlaps between Bloch orbitals at neighboring k- 
points 



M< k < b ) 



ik|«n,k+b) 



(27) 



These overlaps play a central role in the formalism, since 
all contributions to the localization functional can be 
expressed in terms of them. Thus, once the Mi k A b) 
have been calculated, no further interaction with the 
electronic-structure code that calculated the ground state 
wavefunctions is necessary - making the entire Wannier- 
ization procedure a code-independent post-processing 
step There is no unique form for the localization func- 
tional in terms of the overlap elements, as it is possible 



4 Even the case of T-sampling - where the Brillouin zone is sampled 
with a single k-point — is encompassed by the above formulation. 
In this case the neighboring k-points are given by reciprocal lat- 
tice vectors G and the Bloch orbitals differ only by phase factors 
cxp(iG ■ r) from their counterparts at I\ The algebra does be- 
come simpler, thou gh, as will be d iscussed in Sec. |II.F.2| 

5 In particular, see Ferretti et al.\ |2007) for the " 



extension to 



ultrasoft pscudopotcntials and the projector-augmented wave 
method, and |Freimuth et al\ ( |2008[ l; |Kunes et al.\ ( |2010[ ) ; and 
|Posternak et al.\ ( |2002^ for the full-potential linearized aug- 
mented planewave method. 
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to write down many alternative finite-difference expres- 
sions for f„ and (r 2 ) n which agree numerically to leading 
order in the mesh spacing b (first and second order for f„ 
and (r 2 )„ respectively). We give here the expressions of 
Marzari and Vanderbilt (1997), which have the desirable 



property of transforming correctly under gauge transfor- 
mations that shift | On) by a lattice vector. They are 



TV ^ 

k,b 



blm lnM^' b) 



(28) 



(where we use, as outlined in Sec. If. A. 3 the convention 



of Eq. (14 1), and 



k.b 



:,b)|2 



Im In Mj^' b) 



(29) 

The corresponding expressions for the gauge-invariant 
and gauge-dependent parts of the spread functional are 



k,b mn 



(30) 



and 



« = ~E^£i M ^ b) i 2 (31) 
+ ji £ w » £ (- Im lnM - b) - b • f « N 



k,b n 

As mentioned, it is possible to write down alterna- 
tive discretized expressions which agree numerically with 
Eqs. (28)— (31 1 up to the orders indicated in the mesh 
spacing b; at the same time, one needs to be careful in 
realizing that certain quantities, such as the spreads, will 
display slow convergence with respect to the BZ sampling 



(see II. F. 2 for a discussion), or that some exact results 
(e.g., that the sum of the centers of the Wannier func- 
tions is invariant with respect to unitary transformations) 
might acquire some numerical noise. In particular, |Sten-| 
gel and Spaldin ( 2006a ) showed how to modify the above 



expressions in a way that renders the spread functional 
strictly invariant under BZ folding. 



D. Localization procedure 

In order to minimize the localization functional, we 
consider the first-order change of the spread functional 
O arising from an infinitesimal gauge transformation 



TJ<H 



+ dW, 



(k) 



where dW is an infinitesimal 



anti-Hermitian matrix, dW^ — —dW, so that \u n k) 
|«nk) + E m dWml |w mk ) . We use the convention 



dfl \ dfl 



(32) 



(note the reversal of indices) and introduce A and S as 
the superoperators A[B] = {B-B^)/2 and S[B] = (B + 
B^)/2i. Defining 



g( k ' b )=ImlnM^ b )+b-r r , 



o(k,b) _ jiWkjb) TiWk,b> 



T (k,b) 



M (k,b) 



M, 



(k,b) 

k,b) y " 



(33) 



(34) 



(35) 



and referring to |Marzari and Vanderbilt (19971 for the 



details, we arrive at the explicit expression for the gradi- 
ent G^ = dtt/dW^ of the spread functional as 



(36) 



(k) 

ran 



This gradient is used to drive the evolution of the U, 
(and, implicitly, of the |Rn) of Eq. (10)) towards the 
minimum of 57. A simple steepest-descent implementa- 
tion, for example, takes small finite steps in the direction 
opposite to the gradient G until a minimum is reached. 

For details of the minimization strategies and the en- 
forcement of unitarity during the search, the reader is 



referred to Mostofi et al. ( 2008 ) . We should like to point 



out here, however, that most of the operations can be 
performed using inexpensive matrix algebra on small ma- 
trices. The most computationally demanding parts of 
the procedure are typically the calculation of the self- 
consistent Bloch orbitals m the first place, and then 
the computation of a set of overlap matrices 



)u-(0)(k,b) 



(u m \u {0) ) 

\ U mkl"n,k+b/ 



(37) 



that are constructed once and for all from the | ) - Af- 
ter every update of the unitary matrices U^ k ' , the overlap 
matrices are updated with inexpensive matrix algebra 



M (k,b) = ^(kn M (o)(k,b) ^(k+b) 



(38) 



without any need to access the Bloch wavefunctions 
themselves. This not only makes the algorithm compu- 
tationally fast and efficient, but also makes it indepen- 
dent of the basis used to represent the Bloch functions. 
That is, any electronic-structure code package capable 
of providing the set of overlap matrices M( k,b ) can easily 
be interfaced to a common Wannier maximal-localization 
code. 



E. Local minima 

It should be noted that the localization functional can 
display, in addition to the desired global minimum, mul- 
tiple local minima that do not lead to the construction 
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of meaningful localized orbitals. Heuristically, it is also 
found that the WFs corresponding to these local min- 
ima are intrinsically complex, while they are found to be 
real, a part from a single complex phase, at the desired 
global minimum (provided of course that the calculations 
do not include spin-orbit coupling). Such observation in 
itself provides a useful diagnostic tool to weed out unde- 
sired solutions. 

These false minima either correspond to the forma- 
tion of topological defects (e.g., "vortices") in an oth- 
erwise smooth gauge field in discrete k space, or they 
can arise when the branch cuts for the complex loga- 



rithms in Eq. (28) and Eq. (29) are inconsistent, i.e., 



when at any given k-point the contributions from dif- 
ferent b-vectors differ by random amounts of 2tt in the 
logarithm. Since a locally appropriate choice of branch 
cuts can always be performed, this problem is less severe 
than the topological problem. The most straightforward 
way to avoid local minima altogether is to initialize the 
minimization procedure with a gauge choice that is al- 
ready fairly smooth. For this purpose, the projection 
method already described in Sec. |II.B| has been found to 
be extremely effective. Therefore, minimization is usu- 
ally preceded by a projection step, to generate a set of 
analytic Bloch orbitals to be further optimized, as dis- 



Mostofi et al. (2008) 



cussed more fully in Marzari and Vanderbilt (19971 and 



F. The limit of isolated systems or large supercells 

The formulation introduced above can be significantly 
simplified in two important and related cases, which 
merit a separate discussion. The first is the case of open 
boundary conditions: this is the most appropriate choice 
for treating finite, isolated systems (e.g., molecules and 
clusters) using localized basis sets, and is a common ap- 
proach in quantum chemistry. The localization proce- 
dure can then be entirely recast in real space, and corre- 
sponds exactly to determining Foster-Boys localized or- 
bitals. The second is the case of systems that can be 
described using very large periodic supercells. This is 
the most appropriate strategy for non-periodic bulk sys- 
tems, such as amorphous solids or liquids (see Fig. [4] for 
a paradigmatic example), but obviously includes also pe- 
riodic systems with large unit cells. In this approach, the 
Brillouin zone is considered to be sufficiently small that 
integrations over k-vectors can be approximated with a 
single k-point (usually at the T point, i.e., the origin in re- 
ciprocal space). The localization procedure in this second 
case is based on the procedure for periodic boundary con- 
ditions described above, but with some notable simplifi- 
cations. Isolated systems can also be artificially repeated 
and treated using the supercell approach, although care 
may be needed in dealing with the long-range electrostat- 
ics if the isolated entities are charged or have significant 



FIG. 4 (Color online) MLWFs in amorphous Si, either 
around distorted but fourfold coordinated atoms, or in the 
presen ce of a fivefold defect. Adapted from |Fornari et al.\ 

pooil. 



dipole or multipolar character (IDabo et al.\ p2008 ; Makov 



and Payne 1995 ) 



1. Real-space formulation for isolated systems 

For an isolated system, described with open bound- 
ary conditions, all orbitals are localized to begin with, 
and further localization is achieved via unitary trans- 
formations within this set. We adopt a simplified no- 
tation |Rn) — > \w.i) to refer to the localized orbitals of 
the isolated system that will become maximally local- 
ized. We decompose again the localization functional 
n = EJ( r2 h - F ?] into a term Qi = £) a tr [Pr a Qr a ] 
(where P = J2i \ w i)( w i\i Q = 1 — Pi an d 'tr' refers to a 
sum over all the states Wi) that is invariant under unitary 
rotations, and a remainder f2 = J2 a J2i^j \( w i\ r a\wj)\ 2 
that needs to be minimized. Defining the matrices 
Xij = (wi\x\wj), X B ,ij 

similarly for Y and Z), Q can be rewritten as 



XijSij, X' = X - X D (and 



n = tr [X' 2 + Y 



rll 



7/2-1 



(39) 



If X, Y, and Z could be simultaneously diagonalized, 
then would be minimized to zero (leaving only the in- 
variant part). This is straightforward in one dimension, 
but is not generally possible in higher dimensions. The 
general solution to the three-dimensional problem con- 
sists instead in the optimal, but approximate, simulta- 
neous co-diagonalization of the three Hermitian matrices 
X, Y, and Z by a single unitary transformation that 
minimizes the numerical value of the localization func- 
tional. Although a formal solution to this problem is 
missing, implementing a numerical procedure (e.g., by 
steepest-descent or conjugate-gradients minimization) is 
fairly straightforward. It should be noted that the prob- 
lem of simultaneous co-diagonalization arises also in the 



context of multivariate analysis (Flury and Gautschi 



1986) and signal processing (Cardoso and Soulomiac 



1996), and has been recently revisited in relation with 



the present localization approach ( Gygi et al. 2003 ) (see 



also Sec. IIIA in Berghold et al. (2000)). 
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To proceed further, we write 

dQ = 2 tr \X'dX + Y'dY 



Z'dZ] 



(40) 



exploiting the fact that tr [X'Xd] = 0, and similarly 
for Y and Z. We then consider an infinitesimal uni- 
tary transformation \wi) — > \u)i) + dWji\Wj) (where 
dW is antihcrmitian) , from which dX = [X, dW] , and 
similarly for Y, Z. Inserting in Eq. ( |40[ ) and using 
tr[A[B,C}} = tr [C[A,B]] and [X',X] = [X',X U }, we 
obtain dfl = tr [dW G] where the gradient of the spread 
functional G = dfl / dW is given by 



G = 2{[X',X D ] + [Y',Y D ] + [Z',Z D }dw} 



(41) 



The minimization can then be carried out using a proce- 
dure similar to that outlined above for periodic boundary 
conditions. Last, we note that minimizing Q is equiva- 
lent to maximizing tr [X^ + Y§ + Zq], since tr [X'Xd] = 
tr [Y'Yd] = tr [Y'Yd] =0. 



2. T-point formulation for large supercells 

A similar formulation applies in reciprocal space when 
dealing with isolated or very large systems in periodic 
boundary conditions, i.e., whenever it becomes appro- 
priate to sample the wavefunctions only at the T-point 
of the Brillouin zone. For simplicity, we start with the 
case of a simple cubic lattice of spacing L, and define the 
matrices X, y, and Z as 



X„ 



-2irix j L 



Wn) 



(42) 



and its periodic permutations (the extension to supercells 
of arbitrary symmetry has been derived by |Silvestrelli| 
fll999| ))p1 The coordinate x n of the n-th WF center 
(WFC) can then be obtained from 



L 

"2^ 



Im ln(w„|e 



\W„ 



L 

"27T 



Im In X„ 



and similarly for y n and 



(43) 

Equation (43) has been 



shown by Resta (1998) to be the correct definition of 
the expectation value of the position operator for a sys- 
tem with periodic boundary conditions, and had been 
introduced several years ago to deal with the problem 
of determining the average position of a single electronic 
orbital in a periodic supercell (Selloni et al. 1987[ ) (the 



6 We point out that the definition of the X,y,Z matrices for ex- 
tended systems, now common in the literature, is different from 
the one used in the previous subsection for the X, Y, Z matri- 
ces for isolated system. We preserved these two notations for 
the sake of consistency with published work, at the cost of mak- 
ing less evident the close similarities that exist between the two 
minimization algorithms. 



above definition becomes self-evident in the limit where 
w n tends to a Dirac delta function). 

Following the derivation of the previous subsection, 



or of Silvestrelli et al. (19981, it can be shown that the 



maximum-localization criterion is equivalent to maximiz- 
ing the functional 



J 

E 

71=1 



{\x n 



\y n 



| ^"nn | ) 



(44) 



The first term of the gradient d^/dA mn is given by 
[X nm (X* n - X^ m ) - X* nn (X mm - X nn )}, and similarly 
for the second and third terms. Again, once the gradi- 
ent is determined, minimization can be performed using 
a steepest-descent or conjugate-gradients algorithm; as 
always, the computational cost of the localization proce- 
dure is small, given that it involves only small matrices 
of dimension J x J, once the scalar products needed to 
construct the initial X^°\ y^ and have been cal- 
culated, which takes an effort of order J 2 * iVb as is- We 
note that in the limit of a single k point the distinc- 
tion between Bloch orbitals and WFs becomes irrelevant, 
since no Fourier transform from k to R is involved in the 



transformation Eq. (10 1; rather, we want to find the op- 
timal unitary matrix that rotates the ground-state self- 
consistent orbitals into their maximally localized repre- 
sentation, making this problem exactly equivalent to the 
one of isolated systems. Interestingly, it should also be 
mentioned that the local minima alluded to in the pre- 
vious subsection are typically not found when using T 
sampling in large supercells. 

Before concluding, we note that care should be taken 
when comparing the spreads of MLWFs calculated in su- 
percells of different sizes. The Wannier centers and the 
general shape of the MLWFs often converge rapidly as 
the cell size is increased, but even for the ideal case of an 
isolated molecule, the total spread ft often displays much 
slower convergence. This behavior derives from the finite- 
difference representation of the invariant part i of the 
localization functional (essentially, a second derivative); 
while Q i does not enter into the localization procedure, it 
does contribute to the spread, and in fact usually repre- 
sents the largest term. This slow convergence was noted 



by Marzari and Vanderbilt ( 1997 ) when commenting on 



the convergence properties of f2 with respect to the spac- 
ing of the Monkhorst-Pack mesh, and has been studied in 



detail by others ( Stengel and Spaldin 2006a Umari and 



Pasquarello 20031. For isolated systems in a supercell, 



this problem can be avoided simply by using a very large 



L in Eq. ( 43 1 , since in practice the integrals only need to 



be calculated in the small region where the orbitals are 



non-zero (Wu 2004 1. For extended bulk systems, this 
convergence problem can be ameliorated significantly by 
calculating the position operator using real-space inte- 



grals 


Lee 


2006 


Lee et al. 


2005 


Stengel and Spaldin 


2006a 


• 
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G. Exponential localization 



H. Hybrid Wannier functions 



The existence of exponentially localized WFs - i.e., 
WFs whose tails decay exponentially fast - is a famous 
problem in the band theory of solids, with close ties to the 



general properties of the insulating state (Kohn 1964). 
While the asymptotic decay of a Fourier transform can be 
related to the analyticity of the function and its deriva- 



Duffin 


( 


1953 


1 and 



( 1960 ) and references therein) , proofs of exponential lo- 
calization for the Wannier transform were obtained over 
the years only for limited cases, starting with the work of 



Kohn (1959), who considered a one-dimensional crystal 



with inversion symmetry. Other milestones include the 



work of des Cloizeaux ( 1964b ), who established the expo- 



nential localization in ID crystals without inversion sym- 
metry and in the centrosymmetric 3D case, and the sub- 
sequent removal of the requirement of inversion symme- 



try in the latter case by Nenciu (19831. The asymptotic 



behavior of WFs was further clarified by |He and Vande~ 
bilt (2001 1, who found that the exponential decay is mod- 



ulated by a power law. In dimensions d > 1 the problem 
is further complicated by the possible existence of band 
degeneracies, while the proofs of des Cloizeaux and Nen- 
ciu were restricted to isolated bands. The early work on 
composite bands in 3D only established the exponential 



localization of the projection operator P, Eq. (13), not 



of the WFs themselves (des Cloizeaux 1964a). 



The question of exponential decay in 2D and 3D was 



finally settled by Brouder et al. (2007) who showed, as a 



corollary to a theorem by Panati (2007), that a necessary 



and sufficient condition for the existence of exponentially 
localized WFs in 2D and 3D is that the so-called "Chern 
invariants" for the composite set of bands vanish iden- 



tically. Panati (2007) had demonstrated that this con- 
dition ensures the possibility of carrying out the gauge 
transformation of Eq. ^ in such a way that the result- 
ing cell-periodic states \u n k) are analytic functions of k 
across the entire BZQthis in turn implies the exponential 
falloff of the WFs given by Eq. (llOb. 



It is natural to ask whether the MLWFs obtained by 
minimizing the quadratic spread functional Q are also 
exponentially localized. Marzari and Vanderbilt ( 1997 ) 



established this in ID, by simply noting that the MLWF 
construction then reduces to finding the eigenstates of 
the projected position operator PxP, a case for which 



exponential localization had already been proven (Niu 



1991). The more complex problem of exponential local- 



ization of MLWFs for composite bands in 2D and 3D was 



finally proven by Panati and Pisante (2011 ). 



7 Conversely, nonzero Chern numbers pose an obstruction to find- 
ing a globally smooth gauge in k-space. The mathematical defi- 
nition of a Chern number is given in Sec. |V.C.4| 



It is sometimes useful to carry out the Wannier trans- 
form in one spatial dimension only, leaving wavefunc- 
tions that are still delocalized and Bloch-periodic in the 
remaining directions (Sgiarovello et al. 2001). Such or- 
bitals are usually denoted as "hermaphrodite" or "hy- 



brid" WFs. Explicitly, Eq. ( 10 ) is replaced by the hybrid 
WF definition 



Lnkn 



c 

2^ 



2n/c 



k) e 



-ilkj_c 



dk I 



(45) 



where ky is the wavevector in the plane (delocalized di- 
rections) and fcj_, Z, and c are the wavevector, cell index, 
and cell dimension in the direction of localization. The 
ID Wannier construction can be done independently for 
each ky using direct (i.e., non-iterative) methods as de- 
scribed in Sec. IV C 1 of Marzari and Vanderbilt (1997). 



Such a construction has proved useful for a variety of 
purposes, from verifying numerically exponential local- 
ization in 1 dimension, to treating electric polarization 
or applied electric fields along a specific spatial direc- 
tion (iGiustino and Pasquarello 



2005 



Giustino et al. 



[2003} [Murray and Vanderbilt[|2009[|Stengel and Spaldin 
|2006a Wu et al. 2006 ) or for analyzing aspects of topo- 



logical insulators (Coh and Vanderbilt 2009 Soluyanov 
and Vanderbilt 2011a|b ) . Examples will be discussed in 
Sees. [VJL2I and 



VI.A.4I 



I. Entangled bands 

The methods described in the previous sections were 
designed with isolated groups of bands in mind, separated 
from all other bands by finite gaps throughout the entire 
Brillouin zone. However, in many applications the bands 
of interest are not isolated. This can happen, for exam- 
ple, when studying electron transport, which is governed 
by the partially filled bands close to the Fermi level, or 
when dealing with empty bands, such as the four low- 
lying antibonding bands of tetrahedral semiconductors, 
which are attached to higher conduction bands. Another 
case of interest is when a partially filled manifold is to 
be downfoldcd into a basis of WFs - e.g., to construct 
model Hamiltonians. In all these examples the desired 
bands lie within a limited energy range but overlap and 
hybridize with other bands which extend further out in 
energy. We will refer to them as entangled bands. 

The difficulty in treating entangled bands stems from 
the fact that it is unclear exactly which states to choose 
to form J WFs, particularly in those regions of k-space 
where the bands of interest are hybridized with other 
unwanted bands. Before a Wannier localization proce- 
dure can be applied, some prescription is needed for con- 
structing J states per k-point from a linear combination 
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of the states in this larger manifold. Once a suitable J- 
dimensional Bloch manifold has been identified at each 
k, the same procedure described earlier for an isolated 
group of bands can be used to generate localized WFs 
spanning that manifold. 

The problem of computing well localized WFs starting 
from entangled bands is thus broken down into two dis- 
tinct steps, subspace selection and gauge selection. As we 
will see, the same guiding principle can be used for both 
steps, namely, to achieve "smoothness" in k-space. In the 
subspace selection step a J-dimensional Bloch manifold 
which varies smoothly as function of k is constructed. 
In the gauge-selection step that subspace is represented 
using a set of J Bloch functions which are themselves 
smooth functions of k, such that the corresponding WFs 
are well localized. 



1. Subspace selection via projection 



— 



CQ 




FIG. 5 (Color online) Solid black lines: band structure of 
bulk crystalline Si. Blue triangles: band structure for the 
subspace selected by projection onto atomic sp 3 orbitals. Red 
circles: band structure for the subspace selected by projection 
onto atomic sp 3 orbitals and forcing the valence manifold to 
be reproduced exactly using the frozen window indicated. 



The projection technique discussed in Section II. B 



can be easily adapted to produce J smoothly-varying 
Bloch-like states starting from a larger set of Bloch 

bands (|Souza et aL[|2001|). The latter can be chosen, for the symmetry properties of the trial orbitals (|Ku et al. 



often results in "symmetry-adapted" WFs which inherit 



example, as the bands lying within a given energy win- 
dow, or within a specified range of band indices. Their 
number J7k > J is not required to be constant throughout 
the BZ. 

We start from a set of J localized trial orbitals g n (r) 
and project each of them onto the space spanned by the 
chosen eigenstates at each k, 



2002) 



As an example, we plot in Fig.[5]the eight disentangled 
bands obtained by projecting the band structure of sili- 
con, taken within an energy window that coincides with 
the entire energy axis shown, onto eight atomic-like sp 3 
hybrids. The disentangled bands, generated using Wan- 



\4>nk) 



.7k 

E 



\i>mk) (4>mk\9n) 



(46) 



This is identical to Eq. (161, except for the fact that the 
overlap matrix (Ak) mn = (iprnk\9n) has become rectan- 
gular with dimensions J7k x J. We then orthonormalize 



the resulting J orbitals using Eq. (17), to produce a set 



of J smoothly-varying Bloch-like states across the BZ, 



J 

E 



10, 



ikH^k 1/2 )r 



(47) 



As in Eq. (17), (5 k 



(4A) 



-. but 

with rectangular matrices. 

The above procedure achieves simultaneously the two 
goals of subspace selection and gauge selection, although 
neither of them is performed optimally. The gauge se- 
lection can be further refined by minimizing f2 within 
the projected subspace. It is also possible to refine iter- 
atively the subspace selection itself, as will be described 
in the next section. However, for many applications this 
one-shot procedure is perfectly adequate, and in some 
cases it may even be preferable to more sophisticated it- 
erative approaches (see also Sec. III.C ) For example, it 



nier interpolation (Sec. VI.A), are shown as blue trian- 
gles, along with the original first-principles bands (solid 
lines). While the overall agreement is quite good, sig- 
nificant deviations can be seen wherever higher unoccu- 
pied and unwanted states possessing some significant sp 3 
character are admixed into the projected manifold. This 
behavior can be avoided by forcing certain Bloch states 
to be preserved identically in the projected manifold - we 
refer to those as belonging to a frozen "inner" window, 
since this is often the simplest procedure for selecting 
them. The placement and range of this inner window 
will depend on the problem at hand. For example, in or- 
der to describe the low-energy physics for, e.g., transport 
calculations, the inner window would typically include all 
states in a desired range around the Fermi level. 

We show as red circles in Fig. [5] the results obtained 
by forcing the entire valence manifold to be preserved, 
leading to a set of eight projected bands that reproduce 
exactly the four valence bands, and follow quite closely 
the four low-lying conduction bands. For the modifica- 
tions to the projection algorithm required to enforce this 



frozen inner window, we refer to Sec. III.G of Souza et al. 



(2001) 



Projection techniques can work very well, and an ap- 
plication of this approach to graphene is shown in Fig. [6j 
where the it /it* manifold is disentangled with great accu- 
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FIG. 6 (Color online) Solid black lines: band structure of 
graphene. Blue triangles: band structure for the subspace 
selected by projection onto atomic p z orbitals. Red circles: 
band structure for the subspace selected by projection onto 
atomic p z orbitals on each site and sp 2 orbitals on alternate 
sites, and using the frozen window indicated. The lower pan- 
els show the MLWFs obtained from the standard localization 
procedure applied to these two projected manifolds. 



racy by a simple projection onto atomic p z orbitals, or the 
entire occupied manifold together with tt/tt* manifold is 
obtained by projection onto atomic p z and sp 2 orbitals 
(one every other atom, for the case of the sp 2 orbitals - 
albeit bond-centered s orbitals would work equally well). 
Projection methods have been extensively used to 



study strongly-correlated systems 


Anisimov et al. 2005 


Berlijn et al. 


2011 


Ku et al. 


2002 


, in particular to iden- 



tify a "correlated subspace" for LDA+U or DMFT cal- 
culations, as will be discussed in more detail in Sec. |VII| 



It has also been argued (Ku et al. 2010) that projected 



WFs provide a more appropriate basis for the study of 
defects, as the pursuit of better localization in a MLWF 
scheme risks defining the gauge differently for the defect 
WF as compared to the bulk. Instead, a straightforward 
projection approach ensures the similarity between the 
WF in the defect (supercell) and in the pristine (prim- 
itive cell) calculations, and this has been exploited to 
develop a scheme to unfold the band-structure of disor- 
dered supercells into the Brillouin zone of the underlying 
primitive cell, allowing a direct comparison with angu- 



lar resolved photoemission spectroscopy experiments (Ku 



et al. 20101 



2. Subspace selection via optimal smoothness 

The projection onto trial orbitals provides a simple 
and effective way of extracting a smooth Bloch subspace 



starting from a larger set of entangled bands. The rea- 
son for its success is easily understood: the localization 
of the trial orbitals in real space leads to smoothness in 
k-space. In order to further refine the subspace selection 
procedure, it is useful to introduce a precise measure of 
the smoothness in k-space of a manifold of Bloch states. 
The search for an optimally-smooth subspace can then 
be formulated as a minimization problem, similar to the 
search for an optimally-smooth gauge. 

As it turns out, smoothness in k of a Bloch space 
is precisely what is measured by the functional Jli in- 
troduced in Sec. 



II.C.l 



We know from Eq. (19 1 that 
the quadratic spread SI of the WFs spanning a Bloch 
space of dimension J comprises two positive-definite con- 
tributions, one gauge-invariant (Sli), the other gauge- 
dependent (f2). Given such a Bloch space (e.g., an iso- 
lated group of bands, or a group of bands previously dis- 
entangled via projection), we have seen that the opti- 
mally smooth gauge can be found by minimizing Q with 
respect to the unitary mixing of states within that space. 

From this perspective, the gauge-invariance of fii 
means that this quantity is insensitive to the smoothness 
of the individual Bloch states |u nk ) chosen to represent 
the Hilbert space. But considering that fl\ is a part of the 
spread functional, it must describe smoothness in some 
other sense. What Sli manages to capture is the intrinsic 
smoothness of the underlying Hilbert space. This can be 
seen starting from the discretized k-space expression for 



fii, Eq. (30), and noting that it can be written as 



Oi = — 

N 



'k,b 



k,b 



with 



7k,b = Tr[P k Q k+ b], 



(48) 



(49) 



where P k = S n =i l^nkK^nkl * s the gauge-invariant pro- 
jector onto the Bloch subspace at k, Q k = 1 Pu, and 
"Tr" denotes the electronic trace over the full Hilbert 
space. It is now evident that T ki b measures the degree of 
mismatch (or "spillage" ) between the neighboring Bloch 
subspaces at k and k + b, vanishing when they are iden- 
tical, and that fl\ provides a BZ average of the local 
subspace mismatch. 

The optimized subspace selection procedure can now 
be formulated as follows (Souza et al. 2001). A set 



of Jk > J Bloch states is identified at each point on 
a uniform BZ grid, using, for example, an energy win- 
dow. An iterative procedure is then used to extract self- 
consistently across the BZ the J-dimensional subspaces 
having collectively the smallest possible value of SIt (typ- 
ically the minimization starts from an initial guess for 
the target subspaces given, e.g., by projection onto trial 
orbitals). Viewed as function of k, the Bloch subspace 
obtained at the end of the minimization is "optimally 
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smooth" in that it changes as little as possible with k. 
The minimization algorithm can be easily modified in or- 
der to preserve identically the Bloch eigenstates within 
an inner energy window. 

As in the case of the one-shot projection, the outcome 
of this iterative procedure is a set of J Bloch-like states 
at each k which are linear combinations of the initial J7k 
eigenstates. One important difference is that the result- 
ing states are not guaranteed to be individually smooth, 
and the minimization of must therefore be followed 
by a gauge-selection step, which is in every way iden- 
tical to the one described earlier for isolated groups of 
bands. Alternatively, it is possible to combine the two 
steps, and minimize Q — fli + 57 simultaneously with re- 
spect to the choice of Hilbert subspace and the choice 
of gauge (Thygesen et al. 2005a|b ); this should lead to 
the most-localized set of J WFs that can be constructed 
from the initial Jk Bloch states. In all three cases, the 
entire process amounts to a linear transformation taking 
from J7k initial eigenstates to J smooth Bloch-like states, 



\i>nk) 



Jk 

£ 

m— 1 



|VVnk)Vk, ; 



(50) 



In the case of the projection procedure, the explicit ex- 
pression for the J7k x J matrix Vk can be surmised from 



Eqs. (46) and (47) 



Let us compare the one-shot projection and iterative 
procedures for subspace selection, using crystalline cop- 
per as an example. Suppose we want to disentangle the 
five narrow d bands from the wide s band that crosses and 
hybridizes with them, to construct a set of well-localized 
<i-like WFs. The bands that result from projecting onto 
five d-type atomic orbitals are shown as blue triangles in 
Fig. [7] They follow very closely the first-principles bands 
away from the s-d hybridization regions, where the inter- 
polated bands remain narrow. 

The red circles show the results obtained using the 
iterative scheme to extract an optimally-smooth five- 
dimensional manifold. The maximal "global smoothness 
of connection" is achieved by keeping the five d-like states 
and excluding the s-like state. This happens because the 



smoothness criterion embodied by Eqs. (48) and (49) im 



plies that the orbital character is preserved as much as 
possible while traversing the BZ. Inspection of the re- 
sulting MLWFs confirms their atomic d-like character. 
They are also significantly more localized than the ones 
obtained by projection and the corresponding disentan- 
gled bands are somewhat better behaved, displaying less 
spurious oscillations in the hybridization regions. 

In addition, there are cases where the flexibility of the 
minimization algorithm leads to surprising optimal states 
whose symmetries would not have been self-evident in 
advance. One case is shown in Fig. [8] Here we wish to 
construct a minimal Wannier basis for copper, describing 
both the narrow d-like bands and the wide free-electron- 




FIG. 7 (Color online) Solid black lines: band structure of 
bulk crystalline Cu. Blue triangles: band structure for the 
subspace selected by projection onto atomic 3d orbitals. Red 
circles: band structure for the subspace derived from the pre- 
vious one, once the criterion of optimal smoothness has been 
applied. 



like band with which they hybridize. By choosing dif- 
ferent dimensions for the disentangled subspace, it was 
found that the composite set of bands is faithfully repre- 
sented by seven MLWFs, of which five are the standard 
d-like orbitals, and the remaining two are s-like orbitals 
centered at the tetrahedral interstitial sites of the fee 
structure. The latter arise from the constructive interfer- 
ence between sp 3 orbitals that would be part of the ideal 
sp 3 d 5 basis set; in this case, bands up to 20 eV above 
the Fermi energy can be meaningfully described with a 
minimal basis set of seven orbitals, that would have been 
difficult to identify using only educated guesses for the 
projections. 

The concept of a natural dimension for the disentan- 
gled manifold has been explored further by |Thygesen| 
et al. ( 2005a b ) . As illustrated in Fig. |9| they showed 
that by minimizing Q = fli + Q for different choices 
of J, one can find an optimal J such that the result- 
ing "partially-occupied" Wannier functions are most lo- 
calized (provided enough bands are used to capture the 
bonding and anti-bonding combinations of those atomic- 
like WFs). 



3. Iterative minimization of Q% 



The minimization of fli inside an energy window is 



conveniently done using an algebraic algorithm (Souza 
et al. 2001 1. The stationarity condition 5Sli({£t„k}) = 0, 
subject to orthonormality constraints, is equivalent to 
solving the set of eigenvalue equations 



b 



W b fk+b 



I Unit) = A„k | Wnk) • 



(51) 



Clearly these equations, one for each k-point, are 
coupled, so that the problem has to be solved self- 
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FIG. 8 (Color online) Solid black lines: band structure of 
bulk crystalline Cu. Colored lines: band structure for the 
subspace selected by optimal smoothness, and a target di- 
mensionality of 7, giving rise to 5 atom-centered d-like ML- 
WFs and two s-like MLWFs in the tetrahedral interstitials, 
shown below. The color coding represents the projection of 
the disentangled bands onto these MLWFs, smoothly vary- 
ing from red (representing s-like interstitial MLWFs) to blue 
(atom-centered d-like MLWFs). 



the stationary point corresponds to the absolute mini- 
mum of Qj. 



In practice Eq. ( 52 ) is solved in the basis of the origi- 



nal J7k Bloch eigenstates |it„k) inside the energy window. 
Each iteration then amounts to diagonalizing the follow- 
ing J7k x ,7k Hermitian matrix at every k: 



41(k) = U 



b 



(i-i) 



k+b 



(53) 



Since these are small matrices, each step of the iterative 
procedure is computationally inexpensive. 

As a final comment, we mention that in the case of de- 
generacies or quasi degeneracies in the spreads of orbitals 
centered on the same site, the localization algorithm will 
perform rather arbitrary mixing of these (as can be the 
case, e.g., for the d or / electrons of a transition- metal 
ion, or for its ti g or e g groups). A solution to this problem 
is to diagonalize an operator with the desired symmetry 
in the basis of the Wannier functions that have been ob- 
tained (see |Posternak et al. ( 2002 1 for the example of the 
d orbitals in MnO). 



2.85 




0123456789 
Number of unoccupied orbitals {J-N occ ) 



FIG. 9 Plots of average 3 per Wannier function vs. size of the 
Wannier space for several molecules. 3 is defined in Eq. ( |44[ ); 
larger value indicates greater localization. iV occ is the number 
of occupied states, J is the size of the Wannier space, and N 
is the total number of states included in the DFT calculation. 



Adapted from Thygesen et al. (2005a I. 



consistently throughout the Brillouin zone. This can be 
done using an iterative procedure: on the i-th iteration 
go through all the k-points in the grid, and for each of 
them find J orthonormal states defining a sub- 

space whose spillage over the neighboring subspaces from 
the previous iteration is as small as possible. At each step 
the set of equations 



w h P {i - 1] 

w b Mc+b 



-(*) 



X 



(52) 



is solved, and the J eigenvectors with largest eigenvalues 
are selected. That choice ensures that at self-consistency 



4. Localization and local minima 

Empirical evidence and experience suggests that local- 
ized Wannier functions can be readily constructed also 
in the case of an entangled manifold of bands: Even 
for metals, smooth manifolds can be disentangled and 
"wannierized" to give MLWFs. Such disentangled ML- 
WFs, e.g., the p z MLWFs describing the 7r/7r* manifold 
of graphene shown in Fig. |6j are found to be strongly 
localized. While exponential localization has not been 
proven, both numerical evidence and the analogy with 
the isolated composite case suggest this might be the 
case. 

Problems associated with reaching local minima of the 
spread functional, and with obtaining Wannier functions 
that are not real-valued, are more pronounced in the case 
of entangled bands. They are usually alleviated by care- 
ful reconsideration of the energy windows used, in or- 
der to include higher energy states of the appropriate 
symmetry, and/or by using a better initial guess for the 
projections. We infer, therefore, that such problems are 
associated not with the wannierization part of the proce- 
dure, but rather with the initial selection of the smooth 
subspace from the full manifold of entangled bands. 

It is worth noting that the T-point formulation 



(Sec. II.F.2 1 appears to be less affected by these prob- 
lems. In cases where it is not intuitive or obvious what 
the MLWFs should be, therefore, it can often be a fruitful 
strategy to use the T-point formulation to obtain approx- 
imate MLWFs that may then be used to inform the ini- 
tial guess for a subsequent calculation with a full k-point 
mesh. 
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J. Many-body generalizations 

The concept of WFs is closely tied to the framework 
of single-particle Hamiltonians. Only in this case can we 
define J occupied single-particle Bloch functions at each 
wavevector k and treat all J of them on an equal foot- 
ing, allowing for invariance with respect to unitary mix- 
ing among them. Once the two-particle electron-electron 
interaction is formally included in the Hamiltonian, the 
many-body wavefunction cannot be reduced to any sim- 
ple form allowing for the construction of WFs in the usual 
sense. 

One approach is to consider the reduced one-particle 
density matrix 



n(r 



,r') = J **(r,r 2 ...)*(r',r 2 ,...)dr 2 dr 3 .. 



(54) 



for a many-body insulator. Since n(r, r') is periodic in 
the sense of n(r + R, r' + R) = n(r, r'), its eigenvectors - 
the so-called "natural orbitals" - have the form of Bloch 
functions carrying a label n, k. If the insulator is essen- 
tially a correlated version of a band insulator having J 
bands, then at each k there will typically be J occupation 
eigenvalues v n ^ that are close to unity, as well as some 
small ones that correspond to the quantum fluctuations 
into conduction-band states. If one focuses just on the 
subspace of one-particle states spanned by the J valence- 
like natural orbitals, one can use them to construct one- 
particle WFs following the methods described earlier, 



as suggested by Koch and Gocdecker (2001). However, 



while such an approach may provide useful qualitative 
information, it cannot provide the basis for any exact 
theory. For example, the charge density, or expecta- 
tion value of any other one-particle operator, obtained 
by tracing over these WFs will not match their exact 
many-body counterparts. 

A somewhat related approach, adopted by |Hamann] 
and Vanderbilt (2009), is to construct WFs out of the 



quasiparticle states that appear in the GW approxima- 



tion ( Aryasetiawan and Gunnarsson 1998). Such an ap- 



proach will be described more fully in Sec. |VI.A.3[ Here 
again, this approach may give useful physical and chem- 
ical intuition, but the one-electron quasiparticle wave- 
functions do not have the physical interpretation of occu- 
pied states, and charge densities and other ground-state 
properties cannot be computed quantitatively from them. 

Finally, a more fundamentally exact framework for a 
many-body generalization of the WF concept, introduced 



Souza et al. (2000), is to consider a many-body system 



with twisted boundary conditions applied to the many- 
body wavefunction in a supercell. For example, consider 
M electrons in a supercell consisting of M\ x M 2 x M3 
primitive cells, and impose the periodic boundary condi- 
tion 



for j — 1, ...,M, where R is a lattice vector of the su- 
perlattice. One can then construct a single "many-body 
WF" in a manner similar to Eq. Q , but with k — y q and 
iVVik) — > \^ q ) on the right side. The resulting many-body 
WF is a complex function of 3M electron coordinates, 
and as such it is unwieldy to use in practice. However, it 
is closely related to Kohn's theory of the insulating state 



(Kohn, 1964), and in principle it can be used to formu- 



late many-body versions of the theory of electric polar- 
ization and related quantities, as shall be mentioned in 
Sec. rV.A.4l 



III. RELATION TO OTHER LOCALIZED ORBITALS 

A. Alternative localization criteria 

As we have seen, WFs are inherently non-unique and, 
in practice, some strategy is needed to determine the 
gauge. Two possible approaches were already discussed 
in Sec. [TTJ namely, projection and maximal localization. 
The latter approach is conceptually more satisfactory, as 
it does not depend on a particular choice of trial orbitals. 
However, it still does not uniquely solve the problem of 
choosing a gauge, as different localization criteria are pos- 
sible and there is, a priori, no reason to choose one over 
another. 

While MLWFs for extended systems have been gen- 
erated for the most part by minimizing the sum of 
quadratic spreads, Eq. (17), a variety of other localiza- 
tion criteria have been used over the years for molecular 
systems. We will briefly survey and compare some of 
the best known schemes below. What they all have in 
common is that the resulting localized molecular orbitals 
(LMOs) 4>i(r) can be expressed as linear combinations 
of a set of molecular eigenstates ipi(r) (the "canonical" 
MOs), typically chosen to be the occupied ones, 



i (r) = ^C/ ji V J -(r). 



(56) 



The choice of gauge then arises from minimizing or max- 
imizing some appropriate functional of the LMOs with 
respect to the coefficients Uij, under the constraint that 



the transformation ( 56 ) is unitary, which ensures the or- 



thonormality of the resulting LMOs. 

The Foster-Boys criterion (FB) ( |Boys 1960 1966 Fos- 
ter and Boys, 1960a). This is essentially the same as the 



Marzari- Vanderbilt criterion of minimizing the sum of 
the quadratic spreads, except that the sum runs over the 
orbitals in the molecule, rather than in one crystal cell, 



(57) 



i=l 



v I'q(..., Tj + R, ...) = e lq R 5'q(..., Tj, ...) (55) Interestingly, this criterion is equivalent to minimizing 



the "self-extension" of the orbitals (Boys 19661 
E J dridral^^Ol^n-ra) 2 !^^!)! 2 



(58) 



and also to maximizing the sum of the squares of the 
distance between the charge centers of the orbitals 



E K&Mfc 



(59) 



which is closely related to Eq. (44 1 of Sec. II. F. 2 The re- 



lation between Eqs. ( 57 1— ( 59 ) is discussed in Boys (19661. 



The Edmiston-Ruedenberg criterion (ER). Here local- 
ization is achieved by maximizing the Coulomb self- 



interaction of the orbitals (Edmiston and Ruedenbcrg 



1963) 
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(b) 




(c) 



(d) 



FIG. 10 (Color online) (a) One of the two "banana" or r 
orbitals in ethene (b) Centers of the MLWF in ethene (c) r 
MLWF in C0 2 (d) Centers of the MLWF in C0 2 . Atoms 
colors are: hydrogen (white), carbon (yellow), oxygen (red). 
MLWF centers are shown in black. MLWF are computed 
using the scheme of Marzari and Vanderbilt (Marzari and 
1997[) with a large vacuum supercell and a T-point 



Vanderbilt 



sampling of the BZ 



ER. 



E 



d Yl dv 2 |<Mri)| 2 (ri - r 2 )- 1 |0 l (r 2 )| 2 . (60) 



From a computational point of view, the ER approach 
is more demanding (it scales as J 4 versus J 2 for FB), 
but recently more efficient implementations have been 



developed (Subotnik et al. 2007). 



The von Niessen criterion (VN). The goal here is to 



maximize the density overlap of the orbitals ( von Niessen 



1972) 



VN 



J 

E 

i=l 



dndr 2 |^(n)| 2 <J(n - r 2 )|&(r 2 )| 2 . (61) 



The Pipek-Mezey criterion (Pipek and Mezey 1989) 



(PM) . This approach differs from the ones above in that it 
makes reference to some extrinsic objects. The idea is to 
maximize the sum of the squares of the Mulliken atomic 



charges (Mulliken 1955). These are obtained with re- 



spect to a set of atomic orbitals x M centered on atomic 
sites A. We define 



(62) 



where the sum over fj, involves all of the atomic states 
on atom site A, is the density matrix in the atomic 
basis, and — (Xn\Xv) is the overlap operator. The 
functional to be maximized is given by 



J n a 



Q ™ = EE K&i-Wi 



(63) 



= 1 A=l 



This is somewhat similar in spirit to the projection 
scheme discussed in Sec. |II.B| except that it is not a one- 
shot procedure. 



As indicated in Eq. ( 56 ) , all of these schemes amount to 



unitary transformations among the canonical MOs and, 



as such, they give rise to representations of the electronic 
structure of the system that are equivalent to that pro- 
vided by the original set of eigenstates. For the purpose 
of providing chemical intuition, the usefulness of a given 
scheme depends on how well it matches a particular view- 
point of bonding. There have been few studies of the VN 
scheme, but the FB, ER, and PM schemes have been ex- 
tensively compared. In many cases all three approaches 
lead to similar localized orbitals which agree with the 
simple Lewis picture of bonding. A notable exception is 
for systems that exhibit both a and tt bonding. For a 
double bond, both the FB and ER schemes mix the a 
and 7T orbitals, giving rise to two bent "banana bond" 



(or t) orbitals (Pauling 19601 as shown in Fig 10 When 
applied to benzene, both schemes give alternating a and 
r orbitals. Trivially, there are two equivalent sets of these 
orbitals, reminiscent of the two Kekule structures for ben- 
zene. In contrast to the FB and ER schemes, PM pro- 
vides a clear separation of a and tt orbitals. For example, 
in benzene PM gives a framework of six a orbitals and 
three tt orbitals around the ring. 

In some situations the FB and ER orbitals have been 
found to be quite different. This has been observed when 
the bonding in a molecule can be represented as two dis- 
tinct resonant structures. The ER scheme generally gives 
orbitals corresponding to one of the possible structures, 
whilst the FB orbitals correspond a hybrid of the struc- 



tures. An extreme example is C0 2 (Brown et al. 1977). 



In agreement with the 0=C=0 Lewis picture, ER gives 
two lone pairs on each oxygen and two r orbitals between 
each carbon and oxygen. In contrast, the FB scheme 
gives a single lone pair on each oxygen and three highly 
polarized r orbitals between the carbon and each oxygen, 



as shown in Fig. 10 



MLWFs, which may be thought of as the solid-state 
equivalent of FB orbitals, have been widely used to ex- 



19 



amine chemical bonding, as will be discussed in detail in larity measure 



Sec. IV The ER scheme has not been extensively exam- 
ined, an isolated exception being the work of |Miyake and| 
Aryasetiawan ( 2008 ) who proposed a method to maxi- 



mize the Coulomb self-interaction of the orbitals in pe- 



i=l 



\<t>n){<i> n \Ai 



(64) 



riodic systems (see also Sec. VII. B.l). They examined a 
range of bulk transition metals and SrVC>3 and in each 
case found that the resulting WFs were essentially iden- 
tical to MLWFs. 



B. Minimal-basis orbitals 

In the same way as was described for the Marzari 
Vanderbilt scheme of Sec 



with respect to the orthonormal set {</>n}, expressed as 
linear combinations of the original set {ip n }- It is usually 
required that a subset of N < J of the original eigenstates 
are identically preserved ( "frozen in" ) in the disentangled 
subspace, in which case the optimization is performed 
with respect to the remaining J — N states (f> n . J must 
be of sufficient size to capture all of the antibonding char- 
acter of the AOs. 

that 



ILC| the alternative localiza- In later work [t was realized flQian et al.\ |2008[ ) 



tion criteria described above can be applied in a solid- 
state context to isolated groups of bands. One is often 
interested in the more general situation where the bands 
of interest are not isolated. The challenge then is to gen- 
erate a "disentangled" set of J localized WFs starting 
from some larger set of bands. Two procedures for doing 
so were discussed in Sec. II. I namely, projection and it- 



the QOs can be constructed without explicit calculation 
of the eigenstates outside the frozen window. The key 
insight is to realize that the QOs will only have a contri- 
bution from the finite subset of this basis spanned by the 
AOs (see Eq. (67) below). The component of the AOs 
projected onto the N states within the frozen window 
|A|') is given by 



erative minimization of the spread functional f2, not only 
with respect to the choice of gauge, but also with respect 
to the choice of Hilbert space. 

Here we discuss two further procedures which achieve 
the same goal by different means. They have in common 
the fact that the resulting orbitals constitute a minimal 
basis of atomic-like orbitals. 



1. Quasiatomic orbitals 



JV 



i4'> = X>«w»i^>- 



(65) 



The component of \Ai) projected onto the states outside 
the frozen window can hence constructed directly using 
only the AOs and the states within the frozen window, 
as 



I A 



\Ai 



I A 



(66) 



2007| |Lu et al] |2004| |Qian et al.\ |2008[ ) is a projection 



Th e quasiato mic or bitals (QOs ) scheme ( |Chan et al.\ Uging ^ as a basig; the get of {(j>n} which maximize c 



based approach that has the aim of extracting a minimal 
tight-binding basis from an initial first-principles calcu- 
lation. In this regard it is similar in spirit to the Wan- 
nier interpolation techniques discussion in Sec. VI Un- 



like WFs, however, the quasiatomic orbitals form a non- 
orthogonal basis of atom- centered functions, each with 
specified angular character. Their radial part, and hence 
their detailed local shape, depends on the bonding envi- 
ronment. 

As in the Wannier scheme, the first step is to construct 
a suitable J-dimensional subspace {4> n } starting from a 
larger set of J Bloch eigenstates {ipn}- For simplicity 
we consider a T-point only sampling of the BZ and hence 
omit the k index. The goal is to identify a disentangled 
subspace with the desired atomic-orbital character, as 
specified by a given set of J atomic orbitals (AOs) |Aj), 
where i is a composite site and angular-character index. 
One possible strategy would be to employ the one-shot 
projection approach of Sec. |II.I.1| using the AOs as trial 
orbitals. Lu et al. ( 20041) proposed a more optimized pro- 



can be obtained using linear algebra as reported in Qian 

eToTlpoMl ). 

Once a subspace with the correct orbital character has 
been identified, a basis of quasiatomic orbitals can be 
obtained by retaining the portion of the original AOs 
that "lives" on that subspace, 



\Qi) =Y,\<t>n){<l>n\Ai 



(67) 



In general the angular dependence of the resulting QOs 
will no longer that of pure spherical harmonics, but only 
approximately so. In Qian et al. (2008), QOs were ob- 
tained for simple metals and semiconductors. Later ap- 
plications have used the orbitals for the study of quantum 



transport in nano-structures (Qian et al. 20101. 



cedure, based on maximizing the sum-over-square simi- 



8 This is similar in spirit to the Pipek-Mazey procedure, Eq. | |62| 
but applied to subspace selection rather than gauge selection. 
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2. NMTO and Downfolding 

An alternative scheme for obtaining a minimal basis 
representation is the perturbation approach introduced 



by Lowdin ( 1951 1. Here the general strategy is to parti- 
tion a set of orbitals into an "active" set that is intended 
to describe the states of interest, and a "passive" set that 
will be integrated out. Let us write the Hamiltonian for 
the system in a block representation, 



H = 



f H 00 


V 




Voi \ 


V o 


H n J 


> I 


) 



(68) 



where H 00 (i?n) is the projection onto the active (pas- 
sive) subspace, and Vbi is the coupling between the two 
subspaces. An eigenfunction can be written as the sum of 
its projections onto the two subspaces \tp) = + 
This leads to 



(tfoo-e)|Vo> + Voihfc) = 0, 
V 10 \tpo) + (Hu - e)|^i) = 0, 



(69) 
(70) 



where e is the eigenvalue corresponding to Elimina- 
tion of IV'i) gives an effective Hamiltonian for the system 
which acts only on the active subspace: 



H$ (e) = H Q0 - V 01 (Hn - e)" 1 ^ 



10- 



(71) 



This apparent simplification has introduced an energy de- 
pendence into the Hamiltonian. One practical way for- 
ward is to approximate this as an energy-independent 
Hamiltonian Hqq(eo), choosing the reference energy £n 
to be the average energy of the states of interest. 
This approach has been used to construct tight-binding 
Hamiltonians from full electronic structure calculations 



( |Solovyev[|2004 ). 

A form of Lowdin partitioning has been widely used 
in connection with the linear-muffin-tin-orbital (LMTO) 
method, particularly in its most recent formulation, the 



iV th -order mufEn-tin-orbital (NMTO) approach ( 


Ander- 


sen and Saha-Dasgupta 


2000 Zurek et al. 2005 


). In 



this context it is usually referred to as "downfolding" , al- 
though we note that some authors use this term to refer 
to any scheme to produce a minimal basis-set represen- 
tation. 

In studies of complex materials there may be a sig- 
nificant number of MTOs, typically one for each an- 
gular momentum state (s, p, d) on every atomic site. 
One may wish to construct a minimal basis to describe 
states within a particular energy region, e.g., the occu- 
pied states, or those crossing the Fermi level. Let us 
assume we have of basis orbitals (MTOs in this case) 
which we wish to partition into "active" and "passive" 



sets. Using the notation of Eq. (68), Lowdin partitioning 
gives a set of energy-dependent orbitals <po(e,r) for the 
active space according to ( |Zurek et al. 2005 ) 



Taking into account the energy dependence, this reduced 
set of orbitals spans the same space as the original full set 
of orbitals, and can be seen to be the original orbitals of 
the active set dressed by an energy-dependent linear com- 
bination of orbitals from the passive set. In the NMTO 
scheme, the next step is to form an energy-independent 
set of orbitals through an nth-order polynomial fit to the 
energy dependence. 

To give a specific example, accurate calculations on 
tctrahedral semiconductors will require the inclusion of 
d states in the MTO basis, i.e., nine states (s, p, d) per 
site. However, it would be desirable to construct a min- 
imal basis to describe the valence and lower conduction 



states with only four states (s, p) on each site (Lambrecht 



and Andersen 1986 ) . We therefore designate s and p as 
"active" channels and d as "passive". Downfolding will 
result in an MTO with either s or p character on a given 
site, with inclusion of d character on neighboring sites. 
In other words, the tail of the MTO is modified to "fold 
in" the character of the passive orbitals. 



In Fig. 11 the band structure of graphite calculated 
using a full s, p, and d basis on each carbon atom is 



shown in black ( Zurek et al. 2005 1 . The red bands are 



obtained by choosing s, p x , and p y states on every second 
carbon atom as the active channels and downfolding all 
other states. The energy mesh spans the energy range of 
the sp 2 bonding states (shown on the right-hand-side of 



the band-structure plot in Fig. 111. For these bands the 
agreement with the full calculation is perfect to within 
the resolution of the figure. Symmetric orthonormaliza- 
tion of the three NMTOs gives the familiar set of three 
symmetry related cr-bonding orbitals (compare with the 
MLWF of graphene in Fig. [6]). Lechermann et al. (2006) 
have compared MLWFs and NMTO orbitals for a set 
of t 2g states located around the Fermi level SrV03. It 
was found that both schemes gave essentially identical 
orbitals. 



C. Comparative discussion 

At this point it is worth commenting briefly on some 
of the advantages and disadvantages of various choices of 
WFs. We emphasize once again that no choice of WFs, 
whether according to maximal localization or other cri- 
teria, can be regarded as "more correct" than another. 
Because WFs are intrinsically gauge-dependent, it is im- 
possible, even in principle, to determine the WFs experi- 
mentally. By the same token, certain properties obtained 
from the WFs, such as the dipole moments of molecules 



in condensed phases (see Sec. V.B.3 1, must be interpreted 
with caution. The measure of a WF construction proce- 
dure is, instead, its usefulness in connection with theo- 
retical or computational manipulations. 

Where the WFs are to be used as basis functions 



o (e,r) = Mr) - <Mr)(J?u - e) _ Vi, 



(72) for Wannier interpolation (Sec. VI) or other purposes 
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FIG. 11 (Color online) The band structure of graphite calcu- 
lated with a full NMTO spd basis is given in black. The red 
bands have been calculated with an s, p x , and p y orbital on 
every second carbon atom (shown above the band structure). 
Also shown is one of the sp 2 bond orbitals which arise by 
symmetrical orthonormalization of the s, p x , abd p y orbitals. 
The energy meshes used for each calculation are given to the 
right of the band structure. From [Zurek et qT| (|2005). 



(Sec. VII), some variety of maximally localized WFs are 
probably most natural both because the real-space ma- 
trix elements can be restricted to relatively near neigh- 
bors, and because Fourier-transformed quantities become 
relatively smooth in k space. However, in cases in which 
the set of MLWFs does not preserve the space-group sym- 
metry, it may be better to insist on symmetry-preserving 
WFs even at the expense of some delocalization (see also 
the discussion in Sec. II. 1. 1). In this way, properties 



computed from the WFs, such as interpolated bandstruc- 
tures, will have the correct symmetry properties. When 
using WFs to interpret the nature of chemical bonds, as 
in Sec. IV the results may depend to some degree on 
the choice of WF construction method, and the optimal 
choice may in the end be a matter of taste. 



D. Non-orthogonal orbitals and linear scaling 

In recent years there has been much progress in the 
development of practical linear-scaling methods for elec- 
tronic structure calculations, that is, methods in which 
the computational cost of the calculation grows only 
linearly with the size of the system. The fundamen- 
tal principle behind such approaches is the fact that 
electronic structure is inherently local (Heine 1980), or 



'nearsighted' (Kohn 19961. This manifests itself in the 



exponential localization of the WFs in insulators (dis- 



tion properties of the single-particle density matrix 
V 



p(r,r') = (r|P|r') 



(2tt) 3 



BZ 



dk ^/„ k ^nk(r)Ck( r ') ; 



(73) 

where, following Janak (19781, the projection operator 
P of Eq. Q has been generalized to the case of frac- 
tional eigenstate occupancies /„k- The quantity p(r, r') 
has been shown to decay exponentially as exp(— 7|r — r'|) 
in insulators and semiconductors, where the exponent 7 
can be heuristically related to the direct band gap of 



the system ( 


des Cloizeaux 


1964a 


Ismail-Beigi and Arias 


1999 


Taraskin et al. 


2001 


! 1 . It has also been shown to 



with 7 determined by the ratio between the thermal en- 
ergy k^T and the Fermi velocity (Goedecker 1998[ ). 



Exponential localization may seem a surprising result 
given that the Bloch eigenstates extend across the entire 
system. Expressing the density matrix in terms of WFs 
using Eq. (fTol), we find 



p(r,r') = ^ ^ «, iR (r)i^(R' - R>* R ,(r'), (74) 

ij RR' 



Ki 



,(R) = 7^ I dk e-^ R J2i Um ^f^i u{l 



where we have defined the density kerne^®\ 

V 
(2V) 3 

(75) 

which reduces to <5jj<5Ro in the case of a set of fully occu- 
pied bands. We now see that the spatial localization of 
the density matrix is closely linked to that of the Wannier 
functions themselves. This locality is exploited in linear- 
scaling methods by retaining an amount of information in 
the density matrix that scales only linearly with system 
size. 

Many different linear-scaling DFT approaches exist; 



for comprehensive reviews the reader is referred to Galli 
(119961), iGoedeckerl (I1999I), and iBowler and Miyazaki 



(2012). Many of them are based on the variational mini- 
mization of an energy functional expressed either in terms 
of localized Wannier-like orbitals or the density opera- 
tor itself. The common point between these variational 
methods is that the idempotency of the density opera- 
tor or the orthogonality of the Wannier orbitals is not 
imposed explicitly during the minimization procedure. 
Instead, the energy functionals are constructed such that 
these properties are satisfied automatically at the mini- 
mum, which coincides with the true ground state of the 
system. 



9 For metals at zero temperature, the discontinuity in occupancies 
as a function of k results in the well-known algebraic decay of 
the density matrix. 
10 This term was, to the best of our knowledge, first used by 



cussed in Sec. II. G) and, more generally, in the localiza- 



|McWeeny| | Tl960| . 



■MHHBH 



Many of these methods also make use of non- 
orthogonal localized orbitals, referred to as "sup- 
functions" 



(Hernandez and Gillan 19951 



port 

"non-orthogonal generalized Wannier functions" 



or 
(NG- 

WFs) ( |Skylaris et al. 2002), in contrast to canonical 
WFs, which are orthogonal. The density matrix in 
Eq. (74) can be generalized so as to be represented 



in terms of a set of non-orthogonal localized orbitals 
{0aR.( r )} an d a corresponding non-unitary (and, in gen- 
eral, non-square) transformation matrix Af( k \ which 
take the place of {iUiH.(r)} and respectively. Two 

main benefits arise from allowing non-orthogonality. 
First, it is no longer necessary to enforce explicit or- 
thogonality constraints on the orbitals during the en- 



ergy minimization procedure ( Galli and Parrinello 


1992 


Hernandez et al.\ 


1996 


Hierse and Stechel 1994| 


Mauri| 


et al. 


1993 


Ordejon et al. 


1993 


I. Second, a non- 



orthogonal representation can be more localized than an 



essentially equivalent orthogonal one (Anderson 1968 



He and Vanderbilt 2001). In practice, linear-scaling 
methods target large systems, which means that T-point 
only sampling of the BZ is usually sufficient. In this case, 
the separable form for the density matrix simplifies to 

p(r,rO=5> Q (r)i^(r'), (76) 

a/3 

where the density kernel i sp^"| 

=Y}M ] ]Zfn[Mf n . (77) 

n 

Minimization of an appropriate energy functional with 
respect to the degrees of freedom present in the den- 
sity matrix leads to ground-state non-orthogonal orbitals 
that are very similar in appearance to (orthogonal) ML- 
WFs. Fig. [I2| shows, for example, NGWFs on a Ni atom 
in bulk NiO, obtained using the ONETEP linear-scaling 



DFT code (Skylaris et al, 2005). A recent comparison 



of static polarizabilities for molecules, calculated using 



Eq. (89) with both MLWFs and NGWFs, demonstrates 



remarkable agreement between the two (O 'Regan et al. 



2012) 



E. Other local representations 

Over the years, a number of other computational 
schemes have been devised to provide a local analysis 
of the electronic structure in molecules and solids. Here 



11 It is worth noting that the non-orthogonality of the orbitals re- 
sults in an important distinction between covariant and con- 
travariant quantities, as denoted by raised and lowered Greek 
indices l |Artacho and Milans del Bosch] |1991| jO'Regan et al.\ 

[2oTT] i. 



FIG. 12 (Color online) Isosurfaces of the set of nine non- 
orthogonal generalized Wannier functions (NGWFs) on a 
nickel atom in NiO (shown centered on different, symmet- 
rically equivalent, Ni atoms in the lattice). The isosurface is 
set to half of the maximum for the s and p-like NGWFs and 
10 -3 times the maximum for the d-like NGWFs. Adapted 



from O 'Regan et al. (20111 



we briefly mention those most commonly used in solid- 
state studies. The first choice is whether to work with 
the electronic wavefunction or with the charge density. 

One of the earliest and still most widely used wavefunc- 
tion based schemes is the "Mulliken population analysis" 
(Mulliken 1955). This starts from a representation of the 



density operator in an LCAO basis. If an extended basis, 
such as planewaves, has been used, this can be obtained 
after first performing a projection onto a suitable set of 
atomic orbitals (Sanchez-Portal et al. 1995). Using the 



quantity Pa introduced in Eq. (62) the Mulliken charge 
on an atomic site A is given by 



(78) 



The Mulliken scheme also provides a projection into local 
angular-momentum eigenstates and an overlap (or bond) 
population between atom pairs. The major disadvantage 
of the scheme is the fact that the absolute values obtained 
have a marked dependence on the LCAO basis. In fact, 
the results tend to become less meaningful as the basis 
is expanded, as orbitals on one atomic site contribute to 
the wavefunction on neighboring atoms. However, it is 
generally accepted that so long as calculations using the 
same set of local orbitals are compared, trends in the 



values can provide some chemical intuition ( Segall et al. 



1996). An early application was to the study of bonding 



at grain boundaries of Ti02 (Dawson et al. 1996). 



An alternative approach is to work directly with the 



charge density. The scheme of Hirshfeld (1977) attempts 
to partition the charge density by first defining a so-called 
prodensity for the system, typically a superposition of 
free atom charge densities p z (r). The ground-state charge 
density is then partitioned between atoms according to 
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the proportions of the procharge at each point in space. 
This can easily be integrated to give, for example, a total 
charge 



drp(r) 



P*(r) 



(79) 



for each atomic site. Hirschfield charges have recently 
been used to parametrize dispersion corrections to local 
density functionals ( |Tkatchenko and SchefHer 2009). 

Partitioning schemes generally make reference to some 
arbitrary auxiliary system; in the case of Hirschfield 
charges, this is the free-atom charge density, which must 
be obtained within some approximation. In contrast, 
the "atoms in molecules" (AIM) approach developed by 
Baderl (119911) provides a uniquely defined way to partition 



the charge density. It uses the vector field corresponding 
to the gradient of the charge density. In many cases the 
only maxima in the charge density occur at the atomic 
sites. As all field lines must terminate on one of these 
atomic "attractors" , each point in space can be assigned 
to a particular atom according to the basin of attraction 
that is reached by following the density gradient. Atomic 
regions are now separated by zero-flux surfaces S(r s ) de- 
fined by the set of points (r s ) at which 



Vp(r s ) • n(r s ) = 0, 



(80) 



where n(r s ) is the unit normal to S(r s ). Having made 
such a division it is straightforward to obtain values for 
the atomic charges (and also dipoles, quadrupoles, and 
higher moments). The AIM scheme has been widely used 
to analyze bonding in both molecular and solid-state sys- 
tems, as well as to give a localized description of response 



properties such as infra-red absorption intensities ( Matta 



et al. 2007) 



A rather different scheme is the "electron localization 



function" (ELF) introduced by Becke and Edgecombe 



( 1990 ) as a simple measure of electron localization in 
physical systems. Their original definition is based on the 
same-spin pair probability density P(r, r'), i.e., the prob- 
ability to find two like-spin electrons at positions r and 



Savin et al. ( 1992 ) introduced a form for the ELF e(r) 



which can be applied to an independent-particle model: 

1 



e(r) = 



1 + {D/D h f 



1 J 
D = oEN« 



1 ]Vp] 5 



(81) 



(82) 



Dh = ^(3- 2 ) 2/ V /3 ,P = £l^| 2 ' (83) 



i=l 



where the sum is over all occupied orbitals. D repre- 
sents the difference between the non-interacting kinetic 



energy and the kinetic energy of an ideal Bose gas. Dh 
is the kinetic energy of a homogeneous electron gas with 
a density equal to the local density. As defined, e(r) 
is a scalar function which ranges from to 1. Regions 
of large kinetic energy (i.e., electron derealization) have 
ELF values close to zero while larger values correspond 
to paired electrons in a shared covalent bond or in a lone 
pair. In a uniform electron gas of any density, e(r) will 
take the value of 1/2. Early application of the ELF in 
condensed phases provided insight into the nature of the 



bonding at surfaces of Al ( Santis and Resta 2000 ) and 
AI2O3 (Jarvis and Carter 2001), and a large number of 



other applications have appeared since. 



IV. ANALYSIS OF CHEMICAL BONDING 



As discussed in Sec. III.A| there is a long tradition 
in the chemistry literature of using localized molecular 



orbitals ( 


Boys 


19G0 


19GG 


Edmiston and Ruedenberg 


1963 


Foster and Boys 


1960a|b I as an appealing and in- 



tuitive avenue for investigating the nature of chemical 
bonding in molecular systems. The maximally-localized 
Wannier functions (MLWFs) provide the natural gener- 
alization of this concept to the case of extended or solid- 
state systems. Since MLWFs are uniquely defined for the 
case of insulators and semiconductors, they are particu- 
larly suited to discuss hybridization, covalency, and ion- 
icity both in crystalline and disordered systems. In addi- 
tion, in the supercell approximation they can be used to 
describe any disordered bulk, amorphous, or liquid sys- 



tem (Payne et al. 1992), providing a compact descrip- 
tion of electronic states in terms of their Wannier cen- 
ters, their coordination with other atoms, and the spa- 
tial distribution and symmetry of the MLWFs. As such, 
they are often very useful for extracting chemical trends 
and for allowing for statistics to be gathered on the na- 
ture of bonds (e.g covalent bonds vs. lone pairs), be it 
in the presence of structural complexity, as is the case 
of an amorphous solid, or following the intrinsic dynam- 
ics of a liquid or an unfolding chemical reaction. They 
also share the same strengths and weaknesses alluded 
to in section |IILA| whereby different localization crite- 
ria can provide qualitatively different representations of 
chemical bonds. This arbitrariness seems less common 
in extended system, and often some of the most chemi- 
cally meaningful information comes from the statistics of 
bonds as obtained in large-scale simulations, or in long 
first-principles molecular dynamics runs. Finally, local- 
ized orbitals can embody the chemical concept of trans- 
ferable functional groups, and thus be used to construct a 
good approximation for the electronic-structure of com- 
plex systems starting from the orbitals for the different 



fragments (Benoit et al. 2001 Hierse and Stechel 1994 



Lee et al. 2005 ) , as will be discussed in Section VII 
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A. Crystalline solids 

One of the most notable, albeit qualitative, character- 
istics of MLWFs is their ability to highlight the chemical 
signatures of different band manifolds. This was real- 
ized early on, as is apparent from Fig. [2j showing the 
isosurfaces for one of the 4 MLWFs in crystalline silicon 
and gallium arsenide, respectively. These are obtained 
from the closed manifold of four valence bands, yield- 
ing four equivalent MLWFs that map into one another 
under the space-group symmetry operations of the crys- 
tal^ It is clearly apparent that these MLWFs represent 
the intuitive chemical concept of a covalent bond, with 
each MLWF representing the bonding orbital created by 
the constructive interference of two atomic sp 3 orbitals 
centered on neighboring atoms. In addition, it can be 
seen that in gallium arsenide this covalent bond and its 
WFC are shifted towards the more electronegative ar- 
senic atom. This has been explored further by |Abu 



Farsakh and Qteish (2007), to introduce a formal defi- 
nition of electronegativity, or rather of a bond-ionicity 
scale, based on the deviation of WFCs from their geo- 
metrical bond centers. It is worth mentioning that these 
qualitative features of Wannier functions tend to be ro- 
bust, and often independent of the details of the method 
used to obtained them - maximally localized or not. For 
example, similar results are obtained for covalent con- 
ductors whether one makes use of symmetry consider- 



dez et al. 


1997 


Skylaris et al. 


proaches (Stephan et al. 


2000) 



ations (Satpathy and Pawlowska 


1988 


Smirnov et al. 


2002 


Smirnov and Usvyat 2001 


Usvyat et al. 


2004 


), 


finite-support regions in linear scaling methods ( 


Fernan- 



gorithm is applied within Hartree-Fock (Zicovich- Wilson 



et al. 2001 ) rather than density-functional theory. 

Once conduction bands are included via the disen- 
tanglement procedure, results depend on both the tar- 
get dimensions of the disentangled manifold and on the 
states considered in the procedure (e.g., the "outer win- 
dow"). In this case, it becomes riskier to draw conclu- 
sions from the qualitative features of the MLWFs. Still, 
it is easy to see how MLWFs can make the connection 
between atomic constituents and solid-state bands, repre- 
senting a formal derivation of "atoms in solids." That is, 
it can reveal the atomic-like orbitals that conceptually 
lie behind any tight-binding formulation, but that can 
now be obtained directly from first principles according 
to a well-defined procedure. For crystalline silicon, the 



12 It should be noted that the construction procedure does not 
necessarily lead to MLWFs that respect the space-group sym- 
metry. If desired, symmetries can be enforced by imposing co- 
diagonalization of appropriate ope rators (|Posternak et aL||2002 
or by using projection methods ( |Ku et al.\ |2002| |Qian et al. " 
[2008| . 




FIG. 13 (Color online) Four-dimensional manifold (blue cir- 
cles) disentangled from the full band manifold (black lines) for 
Si, together with one of the four antibonding MLWFs that are 
obtained after warmierization of this four-dimensional mani- 
fold. (The other three are equivalent under space-group sym- 
metry operations.) The frozen inner window is also indicated. 



four-dimensional manifold disentangled from the lowest 
part of the conduction bands gives rise to four identical 
antibonding orbitals (see Fig. 13 1 originating from the 
destructive interference of two atomic sp 3 orbitals cen- 
tered on neighboring atoms, to be contrasted with the 
constructive interference shown in Fig. [2] for the valence 
WFs. In addition, an eight-dimensional manifold disen- 
tangled from an energy window including both the va- 
lence bands and the lowest part of the conduction bands 
gives rise to the atomic sp 3 tight-binding orbitals of crys- 
talline silicon (see Fig. 14). These can form the basis of 



the construction of Hamiltonians for model systems (e.g., 
strongly-correlated) or large-scale nanostructures, as will 
be discussed in Chap. |VII| 

These considerations also extend to more complex sys- 
tems. The case of ferroelectric perovskites was studied 



relatively early, as by 


Baranek et al. ( 


20011; Evarestov 


et al. 


(2003 


); and 


Marzari and Vanderbilt 


(19981, thanks 



to the presence of well-separated manifolds of bands 
(King-Smith and Vanderbilt 1994 1. A classic example is 



shown in Fig. |15| showing for BaTiOa the three MLWFs 
per oxygen derived from the nine oxygen 2p bands. While 
in the classical ionic picture of perovskites the B-cation 
(here Ti) is completely ionized in a 4+ state, the covalent 
nature of the bond becomes clearly apparent here, with 
the MLWFs showing a clear hybridization in the form of 
mixed p — d orbitals. Such hybridization is at the origin 
of the ferroelectric instability, as argued by Po sternak 
et aL| ( |1994) . The analysis of the MLWF building blocks 
can also extend to quite different crystal types. For ex- 



ample, Cangiani et al. (2004) discussed the case of Ti02 



polytypes, where bonding MLWFs associated with the O 
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2s /2p valence manifold are seen to be similar in the rutile 
and anatase form, with the third polytype (brookite) an 
average between the two (Posternak et al. 2006). Ap- 
plications to other complex systems can be found, e.g., 
for antiferromagnetic MnO ( Posternak et al. 2002 ) or for 
silver halides ( Evarestov et al. 2004 1 . 



B. Complex and amorphous phases 

Once the electronic ground state has been decom- 
posed into well-localized orbitals, it becomes possible and 
meaningful to study their spatial distribution or the dis- 
tribution of their centers of charge (the WFCs) 



Sil- 



vestrelli et al. (1998) were the first to argue that the 



WFCs can be a powerful tool for understanding bonding 
in low-symmetry cases, representing an insightful and an 
economical mapping of the continuous electronic degrees 
of freedom into a set of classical descriptors, i.e., the cen- 
ter positions and spreads of the WFs. 

The benefits of this approach become apparent when 
studying the properties of disordered systems. In amor- 
phous solids the analysis of the microscopic properties is 
usually based on the coordination number, i.e., on the 
number of atoms lying inside a sphere of a chosen radius 
r c centered on the selected atom (r c is typically inferred 
from the first minimum in the pair correlation function) . 
This purely geometrical analysis is completely blind to 
the actual electronic charge distribution, which ought to 
be important in any description of chemical bonding. An 
analysis of the full charge distribution and bonding in 
terms of the Wannier functions, as for example in Fig. [4] 




FIG. 14 (Color online) Eight-dimensional manifold (blue cir- 
cles) disentangled from the full band manifold (black lines) for 
Si, together with two of the eight atom-centered sp 3 MLWFs 
that are obtained after wannierization of this manifold. (The 
other six are equivalent under space-group symmetry opera- 
tions.) The frozen inner window is also indicated. 



FIG. 15 (Color online) The three maximally-localized Wan- 
nier functions derived from the O 2p bands of BaTiOs, show- 
ing the hybridization with the nominally empty Ti 3d or- 
bitals. The left panel shows the 0[2p z ]-Ti[3rf z2 ] MLWF, the 
central panel adds the 0[2p y ]-Ti[3dj /z ] MLWF, and the right 



panel adds to these the O^p^l-T i^dj; 
from Marzari and Vanderbiltl ( 1998 1 . 



MLWF. Adapted 



for the distorted tetrahedral network of amorphous sili- 
con, would be rather complex, albeit useful to character- 



2001). 



ize the most common defects (Fornari et al. 

Instead, just the knowledge of the positions of the 
WFCs and their spreads can capture most of the chem- 
istry in the system and can identify the defects present. 
In this approach, the WFCs are treated as a second 
species of "classical particles" (representing electrons), 
and the amorphous solid is treated as a statistical assem- 
bly of the two kinds of particles (ions and WFCs). Pair- 
correlation functions can thus be constructed for ions and 
classical electrons, leading to the definition of novel bond- 
ing criteria based on the locations of the WFCs. For the 
case of amorphous silicon, for example, the existence of a 
bond between two ions can be defined by their sharing a 
common WFC within a distance that is smaller than the 
first minimum of the silicon- WFC pair correlation func- 
tion. Following this definition, one can provide a more 
meaningful definition of atomic coordination number, ar- 
gue for the presence (or absence) of bonds in defective 
configurations, and propose specific electronic signatures 



for identifying different defects ( Silvestrelli et al. 1998) 



The ability of Wannier functions to capture the elec- 
tronic structure of complex materials has also been 
demonstrated in the study of boron allotropes. Boron 
is almost unique among the elements in having at least 
four major crystalline phases - all stable or metastable at 
room temperature and with complex unit cells of up to 
320 atoms - together with an amorphous phase. In their 
study of /3-rhombohedral boron, [Ogitsu et al. ( 2009 ) were 
able to identify and study the relation between two-center 
and three-center bonds and boron vacancies, identifying 
the most electron-deficient bonds as the most chemically 
active. Examples are shown in Fig. [16] |Tang and Ismail-] 
Beigi ( 2009 ) were also able to study the evolution of 2D 



boron sheets as they were made more compact (from 
hexagonal to triangular), and showed that the in-plane 



2G 




FIG. 16 (Color online) Charge densities for the MLWFs in /3- 
rhombohedral boron. Red isosurfaces correspond to electron- 
defici ent bonds; b l ue co rrespond to fully occupied bonds. 



From Ogitsu et al. ( 2009 1 





FIG. 17 (Color online) Collapse and amorphization of a Si 
cluster under pressure: increasing to 25 GPa (a) and then to 
35 GPa (b), and then back to 5 GPa (c). Small red "atoms" 
are the Wannier centers. From [Martonak ~et al. (2001). 



bonding pattern of the hexagonal system was preserved, 
with only minor changes in the shape and position of the 
MLWFs. 

Besides its application to the study of disordered net- 



works (Fitzhenry et al. 2003 Lim et al, 2002 Mere- 



galli and Parrinello 2001), the above analysis can also 



be effectively employed to elucidate the chemical and 
electronic properties accompanying structural transfor- 
mations. In work on silicon nanoclusters under pressure 
( |Martonak et al.\ |2000[ |2001| |Molteni et al] |2001| , the 
location of the WFCs was monitored during compressive 
loading (up to 35 GPa) and unloading. Some resulting 
configurations are shown in Fig. [I7j The analysis of the 
"bond angles" formed by two WFCs and their common Si 
atom shows considerable departure from the tetrahcdral 
rule at the transition pressure. The MLWFs also become 
significantly more delocalized at that pressure, hinting at 
a metallization transition similar to the one that occurs 
in Si in going from the diamond to the /3-tin structure. 



C. Defects 



Interestingly, the MLWFs analysis can also point to 
structural defects that do not otherwise exhibit any sig- 
nificant electronic signature 
have predicted - entirely from first-principles 



Goedecker et al. (2002) 



the exis- 
tence of a new fourfold-coordinated defect that is stable 
inside the Si lattice (see Fig. 18). This defect had not 



been considered before, but displays by far the lowest 
formation energy - at the DFT level - among all native 
defects in silicon. Inspection of the relevant "defective" 
MLWFs reveals that their spreads actually remain very 
close to those typical of crystalline silicon, and that the 
WFCs remain equally shared between the atoms in a typ- 
ical covalent arrangement. These considerations suggest 
that the electronic configuration is locally almost indis- 
tinguishable from that of the perfect lattice, making this 
defect difficult to detect with standard electronic probes. 
Moreover, a low activation energy is required for the self- 
annihilation of this defect; this consideration, in combi- 
nation with the "stealth" electronic signature, hints at 
the reason why such a defect could have eluded exper- 
imental discovery despite the fact that Si is one of the 
best studied materials in the history of technology. 

For the case of the silicon vacancy, MLWFs have been 



studied for all the charge states by ( Corsetti and Mostofi 



2011), validating the canonical Watkins model (Watkins 



and Messmer 1974). This work also demonstrated the 



importance of including the occupied defect levels in 
the gap when constructing the relevant WFs, which are 
shown in the first two panels of Fig. [19] For the doubly 
charged split-vacancy configuration, the ionic relaxation 
is such that one of the nearest neighbors of the vacancy 
site moves halfway towards the vacancy, relocating to the 
center of an octahedral cage of silicon atoms. This gives 
rise to six defect WFs, each corresponding to a bond be- 
tween sp 3 d 2 hybrids on the central atom and dangling 
sp 3 orbitals on the neighbors, as shown in the last panel 



of Fig. 19 




FIG. 18 (Color online) The fourfold coordinated defect in Si. 
Si atoms are in green, vacancies in black, and the centers of 
the MLWFs in blue. Adapted from |Goedecker et aZ.| ( |2002| ). 
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FIG. 19 (Color online) Contour-surface plots of the MLWFs 
most strongly associated with a silicon vacancy in bulk sil- 
icon, for different charge states of the vacancy (from left to 
right: neutral unrelaxed, neutral relaxed, and doubly negative 
relaxed). Adapted from Corsetti and Mostofi (20111. 



D. Chemical interpretation 

It should be stressed that a "chemical" interpretation 
of the MLWFs is most appropriate when they are formed 
from a unitary transformation of the occupied subspace. 
Whenever unoccupied states are included, MLWFs are 
more properly understood as forming a minimal tight- 
binding basis, and not necessarily as descriptors of the 
bonding. Nevertheless, these tight-binding states some- 
times conform to our chemical intuition. For example, 
referring back to Fig. |6| we recall that the band struc- 
ture of graphene can be described accurately by disentan- 
gling the partially occupied 7r manifold from the higher 
free-electron parabolic bands and the antibonding sp 2 
bands. One can then construct either a minimal basis of 
one p z MLWF per carbon, if interested only in the ir/n* 
manifold around the Fermi energy, or a slightly larger 
set with an additional MLWF per covalent bond, if in- 
terested in describing both the partially occupied tt/t:* 
and the fully occupied a manifolds. In this latter case, 
the bond-centered MLWFs come from the constructive 
superposition of two sp 2 orbitals pointing towards each 



other from adjacent carbons (Lee et al. 2005 1 . 

On the contrary, as discussed in Sec. |II.I.2| and shown 
in Fig. |8j a very good tight-binding basis for 3d met 



als such as Cu can be constructed ( Souza et al.\ |2001[ ) 
with 5 atom-centered d-like orbitals and two s-like or- 
bitals in the interstitial positions. Rather than reflecting 
a "true" chemical entity, these represent linear combi- 
nations of sp 3 orbitals that interfere constructively at 
the interstitial sites, thus providing the additional vari- 
ational freedom needed to describe the entire occupied 
space. Somewhere in between, it is worth pointing out 
that the atom-centered sp 3 orbitals typical of group-IV 
or III-V semiconductors, that can be obtained in the dia- 
mond/zincblende structure by disentangling the lowest 4 
conduction bands, can have a major lobe pointing either 



to the center of the bond or in the opposite direction ( Lee 



2006 Wahn and Neugebauer 2006). For a given spatial 



cutoff on the tight-binding Hamiltonian constructed from 
these MLWFs, it is found that the former give a quali- 
tatively much better description of the DFT band struc- 



ture than the latter, despite the counter-intuitive result 
that the "off-bond" MLWFs are slightly more localized. 
The reason is that the "on-bond" MLWFs have a sin- 
gle dominant nearest neighbor interaction along a bond, 
whereas for the "off-bond" MLWFs there are a larger 
number of weaker nearest-neighbor interactions that are 



not directed along the bonds (Corsetti 2012). 



E. MLWFs in first-principles molecular dynamics 

The use of MLWFs to characterize electronic bonding 
in complex system has been greatly aided by the imple- 
mentation of efficient and robust algorithms for maximal 
localization of the orbitals in the case of large, and often 
disordered, supercells in which the Brillouin zone can be 
sampled at a single point, usually the zone-center T point 



( Berghold et al. 


2000 


Bernasconi and Madden 


2001 


Sil- 


vestrelli 


1999 


Silvestrelli et al. 


1998 


). Such efforts and 



the implementation in widely-available computer codes 
have given rise to an extensive literature dedicated to 
understanding and monitoring the nature of bonding in 
complex and realistic systems under varying thermody- 
namical conditions or during a chemical reaction. Such 
approaches are particularly useful when combined with 
molecular-dynamics simulations, and most applications 
have taken place within the framework first proposed by 
Car and Parrinello (Car and Parrinello, 19851. In fact, 



specialized algorithms have been developed to perform 
on-the-fly Car-Parrinello molecular dynamics in a Wan- 



nier representation ( 


[ftimie et al. 


2004 


Sharma et al. 


2003 


Wu et al. 


2009 


)■ 



First applications were to systems as diverse as high- 



pressure ice (Bernasconi et al. 1998), doped fullerenes 



vestrelli et al. 



(Billas et al. 



1999), adsorbed organic molecules (Sil 



2000), ionic solids (Bernasconi et al. 2002 



Posternak et al. 2002 1 and the Ziegler-Natta polymeriza- 



tion ( Boero et al. 2000b ) . This latter case is a paradig- 
matic example of the chemical insight that can be gleaned 
by following the WFCs in the course of an first-principles 
simulation. In the Ziegler-Natta reaction we have an in- 
terconversion of a double carbon bond into a single bond, 
and a characteristic agostic interaction between the C-H 
bond and the activated metal center. Both become im- 
mediately visible once the WFCs are monitored, greatly 
aiding the interpretation of the complex chemical path- 
ways. 

Car-Parrinello molecular dynamics is particularly 
suited to the study of liquid systems, and applica- 
tions have been numerous in all areas of physical chem- 



istry. Examples include the work of 


3ako et al.\ 


2002); 


Bernasconi et al. 


2004 


); 


Blumberger et al. (2004); 


Boero 


et al. 


2000a|b); I 


3ucher and Kuyuca 


< ( 


2008); 


Costanzo 


and Delia Valle ( 


2008 


; 


D'Auria et al.\ 


20081 


; van Erp 


and Meijer (2003) 


; Faralli et al. (2006 


;|Heuft and Meijer 


(2005) 


; Ikeda et a 


'. (2005 


); Jungwirth and Tobias 


2002); 
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FIG. 20 (Color online) Snapshots of a rapid water-molecule 
dissociation under high-temperature (1390 K) and high- 
pressure (27 GPa) conditions; one of the MLWFs in the 
proton-donor molecule is highlighted in blue, and one of the 
MLWFs in the proton-acceptor molecule is highlighted in 
green. From Schwegler et al. (2001b I. 



Kirchner and Hutter (20041; Kreitmcir et 



Krekeler et al. (2006); Leung and Rempe (2004); 



and Sprik (2001). 



2003): 



Light- 



stone et al 




2005 1 


>001 


J; Odelius et al. (2004); 


Raugei 


et al. 


(1991 


); 


Raugei and Klein 


(- 


>002); 


Saharay and Bal- 


asubramaniai 


i ( 


200^ 


U; ! 


5alanne et al.\ 


2008) 


; Schwegler 


et al. 


(2001a 


; Sullivan et al. ( 


1999); Suzuki 


(20C 


18); To- 


bias et al. ( 


20011; Todorova et al. 


(2008 


; and 


Vuilleumier 



Water in particular has been studied extensively, both 



at normal conditions (Grossman et al. 2004 Sit et al. 



2007 ) and in high- and low-pressure phases at high tem- 
perature (Boero et al.\ |2000c|d 



2001; 



Romero et al. 



2001] [ Schwegl er et al.\ |2001b| |Silvestrelli and Parrinello 



1999a|b I (a fast dissociation event from one of these sim- 
ulations is shown in Fig. 20 1 . Behavior in the presence of 



et al. 



et al. 



solvated ions (Bako et al] 2002 |Lightstone et al. |2001 



Raugei and Klein[ |2002| |Schwegler et |2001a[ |Tobias 



2001) or a hydrated electron (iBoero 2007 Boero 



2003), or at surfaces and interfaces (Kudin and 



Car 2008 Kuo and Mundy 2004 Kuo et al.\ 2006 2008 



Mundy and Kuo 



2006 



Salvador et all 20031, has also 



been studied. Moreover, MLWFs have been used to cal- 
culate the electronic momentum density that can be mea- 



trato 



2006), and caspases and kinases (Sulpizi et al. 
2003| |2001 1 . Extensive reviews of first-principles quan- 
tum simulations and molecular dynamics, with discus- 
sions of MLWFs in these contexts, have appeared reviews 
by 



Tuckerman (2002); Tuckerman and Martyna 



Dovesi et al. (2005); Kirchner (20071; 



Tse (2002); 



(2000); and 



Vuilleumierl (|2006|), withlMarx and Hutter! (120091) provid 



ing a very comprehensive methodological overview. 

Further applications of first-principles molecular dy- 
namics oriented specifically to extracting information 
about dipolar properties and dielectric responses are dis- 
cussed later in Sec. IV.B.3I 



V. ELECTRIC POLARIZATION AND ORBITAL 
MAGNETIZATION 

First-principles calculations of electric dipoles and 
orbital magnetic moments of molecular systems are 
straightforward. The electric dipole is 



and the orbital moment is 
e \ 



m = — 



2c ^ 

j 



W>j|r x v|V. 



31 > 



(84) 



(85) 



where the sum is over occupied Hamiltonian eigenstates 
r is the position operator, v = (i/ti)[H,r] is the 
velocity operator, and Gaussian units are used. However, 
these formulas cannot easily be generalized to the case of 
crystalline systems, because the Hamiltonian eigenstates 
take the form of Bloch functions Itpnu) that extend over 
all space. The problem is that matrix elements such as 
(ipnk\*\ipnk) an d (ipn\t\r x v|^„k) are ill-defined for such 
extended states (Nenciu |1991 ). 

To deal with this problem, the so-called "modern the- 



work elucidated the relation between the anisotropy of 
the Compton profiles for water and the nature of hydro- 



ory of polarization" (King-Smith and Vanderbilt 1993 



sured in Compton scattering ( |Romero et al] |2000[ ) . This |R es ta[ |1992[ |1994| |Vanderbilt and King-Smith||1993[) was 



developed in the 1990's, and a corresponding "modern 



gen bonding (Romero et al. 2001 ), and led to the sugges- 



tion that the number of hydrogen bonds present can be 



theory of magnetization" in the 2000's (Ceresoli et al. 
Shi eta!] [20071 |Souza and Vanderbnt| |2008| |Thon-| 



directly extracted from the Compton profiles (Sit et al. 



2006 



hauser et al. 



2005 



2007). The population of covalent bond pairs in liq- 



uid silicon and the Compton signature of covalent bond- 



Xiao et al. 2005). Useful reviews of 



these topics have appeared ( 


Rest a 


2000 


2010 


Resta and 


Vanderbilt 


2007 


Vanderbilt and Resta 


2006 


)• 



ing has also recently been studied using MLWFs ( Okada 



et al. 2012) 



These theories can be formulated either in terms of 
Berry phases and curvatures, or equivalently, by working 
in the Wannier representation. The basic idea of the 
latter is to consider a large but finite sample surrounded 
by vacuum and carry out a unitary transformation from 
HIV-1 protease (|Piana et al. 2001), reverse transcrip- the set of delocalized Hamiltonian eigenstates ipj to a 

set of Wannier-like localized molecular orbitals <fij . Then 



Even more complex biochemical systems have been in- 
vestigated, including wet DNA ( Gervasio et al] 2002), 



tase (Sulpizi and Carloni J2000 ), phosphate groups (ATP, 



GTP and ribosomal units) in different environments ( Al 



ber et al. 



19991 |Minehardt et aL||2002"l|Spiegel and Car- 
loni 2003|), drug-DNA complexes (Spiegel and Magis- 



one can use Eq. (84) or Eq. (85), with the tfjj replaced by 



the <f>j , to evaluate the electric or orbital magnetic dipole 
moment per unit volume in the thermodynamic limit. In 
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doing so, care must be taken to consider whether any 
surface contributions survive in this limit. 

In this section, we briefly review the modern theories 
of electric polarization and orbital magnetization and re- 
lated topics. The results given in this section are valid 
for any set of localized WFs; maximally localized ones 
do not play any special role. Nevertheless, the close con- 
nection to the theory of polarization has been one of the 
major factors behind the resurgence of interest in WFs. 
Furthermore, we shall see that the use of MLWFs can 
provide a very useful, if heuristic, local decomposition 
of polar properties in a an extended system. For these 
reasons, it is appropriate to review the subject here. 

A. Wannier functions, electric polarization, and localization 




FIG. 21 Illustration of mapping from physical crystal onto 
equivalent point-charge system with correct dipolar proper- 
ties, (a) True system composed of point ions (+) and charge 
cloud (contours), (b) Mapped system in which charge cloud 
is replaced by quantized electronic charges (— ). In the illus- 
trated model there are two occupied bands, i.e., two Wannier 
functions per cell. 



1. Relation to Berry-phase theory of polarization 

Here we briefly review the connection between the 
Wannier representation and the Berry-phase theory of 



polarization ( 


King-Smith and Vanderbilt 


1994 


Vanderbilt and King-Smith 


1993). 



we have constructed via Eq. ([8|) a set of Bloch-like func- 
tions \ipnk) that are smooth functions of k. Inserting 
these in place of |V>nk) on the right side of Eq. Q, the 
WFs in the home unit cell R=0 are simply 



|0n) 



V 



dk |V> nk ) 



(2^) 3 JBZ 

To find their centers of charge, we note that 



r|0n) = -^3 / rfkHV k e lkr ) 
(2tiT J bz 



|«nk) 



(86) 



(87) 



Performing an integration by parts and applying (0n| on 
the left, the center of charge is given by 



r„ = (OnlrlOn) = 



V 



(2tt) e 



BZ 



dk (M„k|«Vk|Mnk) 



which is a special case of Eq. (23 1. Then, in the home 
unit cell, in addition to the ionic charges +eZ T located 
at positions r T , we can imagine electronic charges — e 
located at positions r n j^] Taking the dipole moment of 
this imaginary cell and dividing by the cell volume, we 
obtain, heuristically 



Z T r T -J2 r * 



(89) 



for the polarization. 



13 In these formulas, the sum over n includes a sum over spin. Al- 
ternatively a factor of 2 can be inserted to account explicitly for 
spin. 



This argument can be put on somewhat firmer ground 
by imagining a large but finite crystallite cut from the 
insulator of interest, surrounded by vacuum. The crys- 
tallite is divided into an "interior" bulk-like region and 
a "skin" whose volume fraction vanishes in the thermo- 
dynamic limit. The dipole moment is computed from 
Eq. (84), but using LMOs <j>j in place of the Hamilto- 
nian eigenfunctions tpj on the right-hand side. Arguing 
that the contribution of the skin to d is negligible in the 
thermodynamic limit and that the interior LMOs become 
bulk WFs, one can construct an argument that arrives 



again at Eq. ( 89 1 



If these arguments still seem sketchy, Eq. ( 89 ) can be 



rigorously justified by noting that its second term 



ol — 
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— V" / dk (u nk \iV k \u nk ) , 
k) 3 „Jbz 



(90) 



is precisely the expression for the electronic contribu- 
tion to the polarization in the Berry-phas e theory ([King- 



Smith and Vanderbilt[ p93l|Resta[[T994l|Vanderbilt and 



King-Smith 1993 ), which was derived by considering the 
flow of charge during an arbitrary adiabatic change of 
the crystalline Hamiltonian. 

The Berry-phase theory can be regarded as provid- 
ing a mapping of the distributed quantum-mechanical 
electronic charge density onto a lattice of negative point 
charges of charge — e, as illustrated in Fig. [21] Then, 
the change of polarization resulting from any physical 
change, such as the displacement of one atomic sublat- 
tice or the application of an electric field, can be related 
in a simple way to the displacements of the Wannier cen- 
ters r„ occurring as a result of this change. 

A well-known feature of the Berry-phase theory is that 
the polarization is only well-defined modulo a quantum 
eR/V, where R is a real-space lattice vector. Such an in- 
determinacy is immediately obvious from Eq. ( 89 ) , since 



the choice of which WFs are assigned to the home unit 
cell (R=0) - or, for that matter, which ions are assigned 
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to it - is arbitrary. Shifting one of these objects by a 
lattice vector R merely changes P by the quantum. Cor- 
respondingly, it can be shown that an arbitrary change 
of gauge can shift individual Wannier centers r„ in ar- 
bitrary ways, except that the sum ^ n r„ is guaranteed 
to remain invariant (modulo a lattice vector). The same 
eH/V describes the quantization of charge transport un- 
der an adiabatic cycle ( |Thouless 1983[ ), and indeed the 
shifts of Wannier charge centers under such a cycle were 
recently proposed as a signature of formal oxidation state 



in crystalline solids (Jiang et al. 2012). 



2. Insulators in finite electric field 

The theory of crystalline insulators in finite electric 
field £ is a subtle one; the electric potential — £ ■ r does 
not obey the conditions of Bloch's theorem, and moreover 
is unbounded from below, so that there is no well-defined 
ground state. In practice one wants to solve for a long- 
lived resonance in which the charge density and other 
properties of the insulator remain periodic, correspond- 
ing to what is meant experimentally by an insulator in a 
finite field. This is done by searching for local minima of 
the electric enthalpy per cell 



F = E KS -V£-P 



(91) 



with respect to both the electronic and the ionic degrees 
of freedom. _Eks is the ordinary Kohn-Sham energy as it 
would be calculated at zero field, and the second term is 
the coupling of the field to the polarization as given in 



Eq. ( 


89 


1 ( 


Nunes and Vanderbilt 


1994 


) or via the 


alent Berry-phase expression ( 


Nunes and Gonze 


Souza et al. 


2002 


Umari and Pasquarello 


2002) 



2001 



This 

approach is now standard for treating periodic insulators 
in finite electric fields in density-functional theory. 



3. Wannier spread and localization in insulators 

We touch briefly here on another interesting connec- 
tion to the theory of polarization. Resta and cowork- 



ers have defined a measure of localization ( 


Resta 2002 


200G 


Resta and Sorella 


1999 


Sgiarovello et al. |2001 



that distinguishes insulators from metals quite generally, 
and have shown that this localization measure reduces, 
in the absence of two-particle interactions or disorder, to 
the invariant part of the spread functional Q\ given in 



Eq. (20). Moreover, Souza et al. (20001 have shown that 



this same quantity characterizes the root-mean-square 
quantum fluctuations of the macroscopic polarization. 
Thus, while the Wannier charge centers are related to the 
mean value of P under quantum fluctuations, their invari- 
ant quadratic spread fli is related to the corresponding 
variance of P. 



4. Many-body generalizations 

In the same spirit as for the many-body WFs discussed 
at the end of Sec. |II. J| it is possible to generalize the for- 
mulation of electric polarization and electron localization 
to the many-body context. One again considers N elec- 
trons in a supercell; for the present discussion we work 
in ID and let the supercell have size L. The many- 
body theory of electric polarization was formulated in 



this context by Ortiz and Martin (1994), and later re- 
formulated by Resta (19981, who introduced a "many- 



body position operator" X = exp(i2-Kx/L) defined in 
terms of the ordinary position operator x = Yli=i^i- 
While (\E r |x|\I f ) is ill-defined in the extended many-body 
ground state |\&), the matrix element | | ^ ) is well- 
defined and can be used to obtain the electric polar- 
ization, up to the usual quantum. These considerations 
were extended to the localization functional, and the re- 
lation between localization and polarization fluctuations, 



by Souza et al. ( 2000 ) . The variation of the many-body 
localization length near an insulator-to-metal transition 
in ID and 2D model systems was studied using quan- 



tum Monte Carlo methods by Hine and Foulkes (2007). 



Finally, the concept of electric enthalpy was generalized 



to the many-body case by Umari et al. (2005), allow- 



ing to calculate for the first time dielectric properties 
with quantum Monte Carlo, and applied to the case of 



the polarizabilities (Umari et al. 2005k and hyperpolariz 



abilities ( Umari and Marzari 2009[ ) of periodic hydrogen 
chains. 



B. Local polar properties and dielectric response 



Is Sec |V.A.l~1 we emphasized the equivalence of the k- 
space Berry-phase expression for the electric polarization, 



Eq. ( 90 ) , and the expression written in terms of the lo- 
cations of the Wannier centers r„, Eq. (89). The lat- 



ter has the advantage of being a real-space expression, 
thereby opening up opportunities for localized descrip- 
tions and decompositions of polar properties and dielec- 
tric responses. We emphasize again that MLWFs have 



no privileged role in Eq. ( 89 ) ; the expression remains 



correct for any WFs that are sufficiently well localized 
that the centers r„ are well defined. Nevertheless, one 
may argue heuristically that MLWFs provide the most 
natural local real-space description of dipolar properties 
in crystals and other condensed phases. 



1. Polar properties and dynamical charges of crystals 

Many dielectric properties of crystalline solids are most 
easily computed directly in the k-space Bloch representa- 
tion. Even before it was understood how to compute the 



polarization P via the Berry-phase theory of Eq. ( 90 1 , it 
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was well known how to compute derivatives of P using 



linear-response methods (Baroni et al. 2001 de Gironcoli 



et al] |1989[ |Resta| |1992[ ) . Useful derivatives include the 
electric susceptibility Xij = dPi/dEj and the Born (or 
dynamical) effective charges Z* T ^ = VdPi /dR T j , where i 
and j are Cartesian labels and R T j is the displacement of 
sublattice r in direction j. With the development of the 
Berry-phase theory, it also became possible to compute 
effective charges by finite differences. Similarly, with the 



electric-enthalpy approach of Eq. (91 1 it became possible 



to compute electric susceptibilities by finite differences as 



well (Souza et al. 2002 Umari and Pasquarello 2002) 



The Wannier representation provides an alternative 
method for computing such dielectric quantities by fi- 
nite differences. One computes the derivatives dr n ^/dEj 
or dr n ^/ dR T j of the Wannier centers by finite differences, 
then sums these to get the desired Xij or ^trj- An exam_ 
ple of such a calculation for Z* in GaAs was presented 



already in Sec. VII of Marzari and Vanderbilt 



and an application of the Wannier approach of 



1997), 



Nunes| 



and Vanderbilt (1994) in the density-functional context 



was used to compute x by finite differences for Si and 
GaAs (Fernandez et al. 1998). Dynamical charges were 



computed for several Ti0 2 phases by |Cangiani et al.\ 
(2004) and Posternak et al. (20061, and, as mentioned 



in Sec. |IV.A observed differences between polymorphs 
were correlated with changes in the chemical nature of 
the WFs associated with OTi3 structural units. Piezo- 
electric coefficients, which are derivatives of P with re- 
spect to strain, have also been carried out in the Wannier 



representation for ZnO and BeO by Noel et al. (2002) 



Some of the most extensive applications of this kind 
have been to ferroelectric perovskites, for which the 
dynamical charges have been computed in density- 
functional and/or Hartrce-Fock contexts for BaTi0 3 , 
KNb0 3 , SrTi0 3 , and PbTi0 3 dBaranek et al. " 



Evarestov et ail 12003 Marzari and Vanderbilt 



2001 



1998 



Usvyat et al. 2004). In these materials, partially covalent 



bonding associated with hybridization between O 2p and 
Ti 3d states plays a crucial role in stabilizing the ferro- 
electric state and generating anomalous dynamical effec- 



tive charges ( Posternak et al.\ 1994 Zhong et al. 1994). 



Recall that the dynamical, or Born, effective charge Z* is 
defined as Z* T j = VdPi/dR r j and carries units of charge. 
Naively, one might expect values around +4e for Ti ions 
and — 2e for oxygen ions in BaTi0 3 based on nominal ox- 
idation states, but instead one finds "anomalous" values 



that are much larger. For example, Zhong et al. (1994) 
reported values of +7.2e for Ti displacements, and 
for O displacements along the Ti-O-Ti chains. 



(a) 



(b) 




FIG. 22 (Color online) Amplitude isosurface plots of the 
maximally-localized 0[2p z ]-Ti[3G? z 2] Wannier functions in 
BaTiOs. O is at center, surrounded by a plaquette of four Ba 
atoms (green); Ti atoms (yellow, almost hidden) are above 
and below, (a) Centrosymmetric structure, (b) Ferroelectric 
structure in which central O is displace upward relative to 
neighboring Ti atoms. 



the changes in the MLWFs induced by the atomic dis- 
placements. Fig. |22"t» shows an Q[2p z ]-Ti[3d z 2] MLWF 



in centrosymmetric BaTi0 3 (Marzari and Vanderbilt 



1998). The hybridization to Ti 3d z 2 states appears in 



the form of the "donuts" surrounding the neighboring 
Ti atoms. When the O atom moves upward relative to 
the geometric center of the two neighboring Ti atoms as 



shown in Fig. 22 b), as it does in ferroelectrically dis- 
torted BaTi0 3 , the hybridization strengthens for the up- 
per O-Ti bond and weakens for the lower one, endowing 
the WF with more Ti 3c? character on the top than on the 
bottom. As a result, the center of charge of the WF shifts 
upward, and since electrons carry negative charge, this 
results in a negative anomalous contribution to the Z* 
of the oxygen atom. The figure illustrates this process 
for (7-oriented oxygen WFs, but a similar effect occurs 
for the 7r-oriented oxygen WFs, and the total anomalous 
dynamical charge can be accounted for quantitatively on 
the basis of the distortion-induced changes of each kind 



of WF in the crystal (Marzari and Vanderbilt 1998) 



The above illustrates the utility of the MLWFs in pro- 
viding a local description of dielectric and polar responses 
in crystals. This strategy can be carried further in many 
ways. For example, it is possible to decompose the Z* 
value for a given atom in a crystal into contributions 
coming from various different neighboring WFs, as was 
done for GaAs in Sec. VII of Marzari and Van derbiltl 



-5.7e 



( 1997 ) and for BaTi0 3 by Marzari and Vanderbilt ( 1998 ) . 



Some chemical intuition is already gained by carrying out 
a band-by-band decomposition of the Z* contributions 
This behavior can be understood (Posternak et al. (Ghosez and Gonze 2000 Ghosez et al. 19951, but the 



1994 Zhong et al. 1994 ) as arising from hybridization be- 



tween neighboring O 2p and Ti 3c? orbitals that dominate 
the valence and conduction bands, respectively. This 
hybridization, and the manner in which it contributes 
to an anomalous Z*, can be visualized by inspecting 



WF analysis allows a further spatial decomposition into 
individual WF contributions within a band. A deeper 
analysis that also involves the decomposition of the WFs 
into atomic orbitals has been shown to provide further in- 



sight into the anomalous Z* values in perovskites (Bhat 
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FIG. 23 (Color online) MLWFs for /3-PVDF polymer chain, 
(a) MLWF charge centers, indicated by diamonds, (b)-(d) 
MLWFs localized on C-C, C-F, and C-H bonds, respectively. 
From 



Nakhmanson et al. (2005) 



„ 0.8- 



0.6- 



I 0.4- 



0.2- 



0- 




FIG. 24 (Color online) Dispersion of WF center positions 
along z as a function of (k x ,k y ) for a superlattice composed 
of alternating layers of SrTiOs (sublayers 1A and IB) and 
BaTi0 3 (sublayers 2A and 2B). From|Wu eFaTI (|2006l. 



tacharjee and Waghmare, 2010) 



Some insightful studies of the polar properties of poly- 
mer systems in terms of MLWFs have also been carried 



out. Figure [23j for example, shows the WF centers and 
characters for the /3 conformation of polyvinylidene flu- 



oride (/3-PVDF) ( |Nakhmanson et al] |2005[ ), one of the 
more promising ferroelectric polymer systems. An in- 
spection of WF centers has also been invoked to explain 
the polar properties of so-call "push-pull" polymers by 



interfaces and related systems. Later applications to per- 
ovskite oxide films and superlattices have been fairly ex- 
tensive. The essential observation is that, when studying 
a system that is layered or stacked along the z direction, 
one can still work with Bloch functions labeled by k x and 
k y while carrying out a Wannier construction or analysis 
only along z. Since the extraction of Wannier centers in 



Kudin et al.\ ([2007) and of H 2 ice by ILu et al.\ p08|. ID is rather trivial, even in the multiband case flBhat 



Finally, we note an interesting recent study in which 
changes in polarization induced by corrugations in BN 



sheets were analyzed in terms of WFs ( |Naumov et al. 
20091). 



2. Local dielectric response in layered systems 

In a similar way, the theoretical study of dielectric 
properties of ultrathin films and superlattices can also 
be enriched by a local analysis. Two approaches have 
been introduced in the literature. In one, the local x-y- 
averaged electric field £ z (z) is calculated along the stack- 
ing direction z, and then the local dielectric permittiv- 
ity profile e(z) — £ z {z)/D z or inverse permittivity pro- 
file e~ 1 (z) = D z /£ z (z) is plotted, where D z is the x-y- 
averaged electric displacement field (constant along z in 
the absence of free charge) determined via a Berry-phase 
calculation of P z or by inspection of £ z in a vacuum re- 
gion. Such an approach has been applied to study di- 
electric materials such as S1O2 and Hf02 interfaced to 



Si( 


Giustino and Pasquarello 


2005 


Giustino et al. |2003| 


Shi and Ramprasad 


2006 


2007 1 and perovskite films and 


superlattices ( Stengel and Spaldin 


2006a|b 


)• 



The second approach is to use a Wannier analysis to as- 
sign a dipole moment to each layer. This approach, based 
on the concept of hybrid WFs discussed in Sec. |II.H| 



was pioneered by Giustino and Pasquarello ( 2005[ ) and 



tacharjee and Waghmare} 2005| Marzari and Vandcrbilt 



1997 Sgiarovello et al. 2001), it is not difficult to con- 
struct a "Wannier center band structure" z(k x ,k y ), and 
use the planar-averaged values to assign dipole moments 
to layers. This approach was demonstrated by |Wu et al. 
(2006), as shown in Fig. 24 and has since been used to 



study perovskite superlattices and artificial nanostruc- 



tures (Murray and Vanderbilt 2009 Wu et al. 20081. 



In a related development, Stengel and Spaldin (2006a 



introduced a Wannier-based method for computing po- 
larizations along z and studying electric fields along z 
that work even in the case that the stacking includes 
metallic layers, as long as the system is effectively insu- 
lating along z (Stengel and Spaldin 2007). This allows 



for first-principles calculations of the nonlinear dielectric 
response of ultrathin capacitors and related structures 
under finite bias, providing an insightful avenue to the 
study of finite-size and dead-layer effects in such systems 
( [Stengel et aT\ |2009a|b|c I . 



Giustino et al. (2003) and used by them to study Si/Si02 



3. Condensed molecular phases and solvation 

Wannier-function methods have also played a promi- 
nent role in the analysis of polar and dielectric proper- 
ties of dipolar liquids, mainly H 2 and other H-bonded 
liquids. While the dipole moment of an isolated H 2 
molecule is obviously well defined, a corresponding defi- 
nition is not easy to arrive at in the liquid phase where 
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molecules are in close contact with each other. An influ- 
ential development was the proposal made by |SilvestreTE 



and Parrinello ( 1999a|b ) that the dipole moment of a 
solvated molecule could be defined in terms of positive 
charges on ionic cores and negative charges located at the 
centers of the MLWFs. Using this definition, these au- 
thors found that the water molecule dipole is somewhat 
broadly distributed, but has an average value of about 
3.0 D, about 60% higher than in the gas phase. These 
features were shown to be in conflict with the behavior 
of widely-used empirical models. 

Of course, such a definition in terms of the dipole of the 
molecular WF-center configuration remains heuristic at 
some level. For example, this local measure of the dipole 
does not appear to be experimentally measurable even 
in principle, and clearly the use of one of the alternative 
measures of maximal localization discussed in Sec. IIII.AI 
would give rise to slightly different values. Nevertheless, 
the approach has been widely adopted. For example, 



subsequent work elaborated ( Dyer and Cummings 2006 



Sagui et al. 2004 ) and extended this analysis to water 



in supercritical conditions (Boero et al. 2000c|d l, con 



2005), and with solvated ions present (Scipioni et al. 



fined geometries ( Coudert et al. 2006 Dellago and Naor 



2009), and compared the results obtained with different 



exchange-correlation functionals ( Todorova et al. 2006 1 . 

It should be noted that the decomposition into Wan- 
nier dipoles is closer to the decomposition of the charge 
density into static (Szigeti) charges than to a decom- 
position into dynamical (Born) charges. The first one 
corresponds to a spatial decomposition of the total elec- 
tronic charge density, while the second is connected with 
the force that appears on an atom in response to an ap- 
plied electric field. As a counterpoint to the WF-based 



definition, therefore, Pasquarello and Resta (2003) have 



argued the a definition based on these forces provides a 
more fundamental basis for inspecting molecular dipoles 
in liquids. In particular, they define a second-rank ten- 
sor as the derivative of the torque on the molecule with 
respect to an applied electric field, and finding that this 
is typically dominated by its antisymmetric part, they 
identify the latter (rewritten as an axial vector) as the 
molecular dipole. Surprisingly, they find that the mag- 
nitude of this dipole vector is typically only about 2.1 D 
in liquid water, much less than the value obtained from 
the WF analysis. 

Clearly the WF-based and force-based approaches to 
defining molecular dipoles provide complementary per- 
spectives, and a more complete reconciliation of these 
viewpoints is the subject of ongoing work. 

Finally, we note that there is an extensive literature 
in which Car-Parrinello molecular-dynamics simulations 
are carried out for H 2 and other liquids, as already 



( Kirchner 


2007 


Tse[|2002; 


Tuckerman 


2002 


Tuckerman 


and Martyna 


2000 


). Using such approaches, it is possi- 



ble to compute the dynamical dipole-dipole correlations 
of polar liquids and compare the results with experimen- 
tal infrared absorption spectra. While it is possible to 
extract the needed information from the time-time corre- 
lation function of the total polarization P(t) of the entire 
supercell as calculated using the Berry-phase approach, 
methods which follow the time-evolution of local dipoles, 
as defined via WF-based methods, provide additional in- 



sight and efficiency ( Bernasconi et al. 


1998 


Chen et al. 


2008 McGrath et al. 2007| Pasquarel 


o and 


ilesta[ 


2003 


Sharma et al. 2005 2007 ) and can easily be extended to 


other kinds of molecular systems (Gaigeot et al.\ 


2007 



Gaigeot and Sprikj [20031 |Gaigeot et al.\ |2005| |Pagliai 



et al. 2008 ) . The applicability of this kind of approach 
has benefited greatly from the development of methods 
for computing WFs and their centers "on the fly" during 



Car-Parrinello molecular dynamics simulations (Iftimic 



aL 2004 Sharma et al. 2003 Wu et al. 2009). 



C. Magnetism and orbital currents 

1. Magnetic insulators 

Just as an analysis in terms of WFs can help clarify 
the chemical nature of the occupied states in an ordi- 
nary insulator, they can also help describe the orbital 
and magnetic ordering in a magnetic insulator. 

In the magnetic case, the maximal localization pro- 
ceeds in the same way as outlined in Sec. [XXJ, with trivial 
extensions needed to handle the magnetic degrees of free- 
dom. In the case of the local (or gradient-corrected) spin- 
density approximation, in which spin-up and spin-down 
electrons are treated independently, one simply carries 
out the maximal localization procedure independently 
for each manifold. In the case of a spinor calculation 
in the presence of spin-orbit interaction, one instead im- 
plements the formalism of Sec. [IT] treating all wavefunc- 
tions as spinors. For example, each matrix element on 



the right-hand side of Eq. (27) is computed as an inner 
product between spinors, and the dimension of the re- 
sulting matrix is the number of occupied spin bands in 
the insulator. 

Several examples of such an analysis have appeared 
in the literature. For example, applications to simple 



antiferromagnets such as MnO ( jPosternak et al. 2002 \ , 
novel insulating ferromagnets and antiferromagnets (Ku 



et al. 2002 2003), and complex magnetic ordering in 



rare-earth manganates (Picozzi et al. 2008 Yamauchi 



et al. 2008 ) have proven to be illuminating. 



2. Orbital magnetization and NMR 

In a ferromagnetic (or ferrimagnetic) material, the to- 
tal magnetization has two components. One arises from 
electron spin and is proportional to the excess population 
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of spin-up over spin-down electrons; a second corresponds 
to circulating orbital currents. The spin contribution is 
typically dominant over the orbital one (e.g., by a factor 
of 10 or more in simple ferromagnets such as Fe, Ni and 



Co (Ceresoli et al. 2010a)), but the orbital component 
is also of interest, especially in unusual cases in which it 
can dominate, or in the context of experimental probes, 
such as the anomalous Hall conductivity, that depend on 
orbital effects. Note that inclusion of the spin-orbit inter- 
action is essential for any description of orbital magnetic 
effects. 

Naively one might imagine computing the orbital mag- 
netization M or b as the thermodynamic limit of Eq. ( 85 ) 
per unit volume for a large crystallite. However, as we 
discussed at the beginning of Sec. [Vj Bloch matrix ele- 
ments of the position operator r and the circulation oper- 
ator r x v are ill-defined. Therefore, such an approach is 
not suitable. Unlike for the case of electric polarization, 
however, there is a simple and fairly accurate approxi- 
mation that has long been used to compute M or b: one 
divides space into muffin-tin spheres and interstitial re- 
gions, computes the orbital circulation inside each sphere 
as a spatial integral of r x J, and sums these contribu- 
tions. Since most magnetic moments are fairly local, such 
an approach is generally expected to be reasonably accu- 
rate. 

Nevertheless, it is clearly of interest to have available 
an exact expression for M or b that can be used to test 
the approximate muffin-tin approach and to treat cases 
in which itinerant contributions are not small. The solu- 
tion to this problem has been developed recently, leading 
to a closed-form expression for M or b as a bulk Brillouin- 
zone integral. Derivations of this formula via a semi- 



classical (Xiao et al. 2005) or long- wave quantum (Shi 



et al. 2007) approach are possible, but here we empha- 



size the derivation carried out in the Wannier represen- 



tation ( 


Ceresoli et al. 


2006 


Souza and Vanderbilt 


2008 


Thonhauser et al. 


2005 


). For this purpose, we restrict 



our interest to insulating ferromagnets. For the case of 
electric polarization, the solution to the problem of r ma 



where Aq is the unit cell area, since in the interior the 
LMOs (j>j are really just bulk WFs. This time, however, 
the skin contribution does not vanish. The problem is 
that (4>j | v 1 j ) is nonzero for LMOs in the skin region, and 
the pattern of these velocity vectors is such as to describe 
a current circulating around the boundary of the sample, 
giving a second "itinerant circulation" contribution that 
can, after some manipulations, be written in terms of 
bulk WFs as 



)]. (93) 



Mic = -^-z 2JJ2x<R|y|0) - RyWO 

R 



When these contributions are converted back to the 
Bloch representation and added together, one finally ob- 
tains 



M, 



orb 



2hc 



Im 



d 2 k 

W) 2 



(<9 k Uk|x (H u + £ k )|<9 k u k ) 



which is the desired k-space bulk expression for the or- 
bital magnetization ( Thonhauser et al. 2005 ) ^] The 
corresponding argument for multiple occupied bands in 



three dimensions follows similar lines (Ceresoli et al. 



2006 Souza and Vanderbilt 2008), and the resulting for- 



mula has recently been implemented in the context of 



pscudopotential plane- wave calculations (Ceresoli et al. 



2010a). Interestingly, it was found that the interstitial 
contribution - defined as the difference between the to- 



tal orbital magnetization, Eq. (94), and the muffin-tin 
result - is not always negligible. In bec Fe, for example, 
it amounts to more than 30% of the spontaneous orbital 
magnetization, and its inclusion improves the agreement 
with gyromagnetic measurements. 

The ability to compute the orbital magnetization is 
also of use in obtaining the magnetic shielding of nu- 
clei. This is responsible for the chemical shift effect 
observed in nuclear magnetic resonance (NMR) exper- 
iments. A first principles theory for magnetic shielding 
in solids was established by examining the perturbative 
response to a periodic magnetic field in the long wave- 



trix elements was sketched in Sec. 17X11 and a heuristic ^ h Umit 1 Mauri et al \ l 1996a l l Pickard and Mauri 



derivation of Eq. ( 89 ) was given in the paragraph follow- 
ing that equation. A similar analysis was given in the 
above references for the case of orbital magnetization, as 
follows. 

Briefly, one again considers a large but finite crystallite 
cut from the insulator of interest, divides it into "inte- 
rior" and "skin" regions, and transforms from extended 
eigenstates to LMOs cj)j. For simplicity we consider the 
case of a two-dimensional insulator with a single occupied 
band. The interior gives a rather intuitive "local circula- 
tion" (LC) contribution to the orbital magnetization of 
the form 



2001). An alternative perturbative approach used a WF 



representation of the electronic structure together with a 



periodic pos ition operator (|Sebastiani| 



et al. 2002 Sebastiani and Parrinello 



20031 [Sebastiani 



2001). However, 



magnetic shieldings can also be computed using a "con- 



verse" approach in which one uses Eq. (94) to compute 



the orbital magnetization induced by a fictitious point 



magnetic dipole on the nucleus of interest ( Ceresoli et al. 



Ah 



LC 



2,4 c 



(0|r x v|0) 



(92) 



14 In the case of metals Eq. \9A\ must be modified by adding a 
—Ijx term inside the parenthesis, with /i the chemical potential 
( |Ceresoli et al.\ |2006[ |Xiao et a7T]|2005| |. Furthermore, the inte- 
gration is now restricted to the occupied portions of the Brillouin 
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2010b Thonhauser et al. 2009). The advantage of such insulator, and Thonhauser and Vanderbilt (20061 have 



approach is that it does not require linear-response the- 
ory, and so it is amenable to large-scale calculations or 
complex exchange-correlation functionals (e.g., including 
Hubbard U corrections, or Hartree-Fock exchange), al- 
beit at the cost of typically one self-consistent iteration 
for every nucleus considered. Such converse approach has 
then been extended also to the calculation of the EPR 



g-tensor by Ceresoli et al. (2010a 



3. Berry connection and curvature 

Some of the concepts touched on in the previous sec- 
tion can be expressed in terms of the k-space Berry con- 
nection 



explored the way in which the usual Wannier construc- 
tion breaks down for model systems. 

Second, depending on how their Bloch functions wrap 
the Brillouin zone, nonmagnetic (T-invariant) insulators 
can be sorted into two classes denoted as "i^-even" and 
"Z 2 -odd" (after the name Z 2 of the group {0, 1} un- 
der addition modulo 2). Most (i.e., "normal") insulators 
are i^-even, but strong spin-orbit effects can lead to the 
2 2 -odd state, for which the surface-state dispersions are 
topologically required to display characteristic features 
that are amenable to experimental verification. Several 
materials realizations of Z^-odd insulators have now been 



confirmed both theoretically and experimentally (Hasan 



and Kane 2010 Hasan and Moore 2011). 



A„ k = (Mnk|«V k |u„k) 

and Berry curvature 

•F nk = x A n k 



(95) 



(96) 



of band n. In particular, the contribution of this band 



to the electric polarization of Eq. (l90l, and to the second vious subsection ( |Essin et al.\ |2009| |Qi et al.\ |2008[ ) 



In a related development, the orbital magnetoelectric 
coefficient aty = 9Af or bj /d£i was found to contain an 
isotropic contribution having a topological character (the 
"axion" contribution, corresponding to an £ ■ B term in 
the effective Lagrangian). This term can be written as a 
Brillouin-zone integral of the Chern-Simons 3-form, de- 
fined in terms of multiband generalizations of the Berry 
connection Ak and curvature introduced in the pre- 

The 



term in the orbital magnetization expression of Eq. (94), 



are proportional to the Brillouin-zone integrals of A, ik 
and -Enk-^nk, respectively. These quantities will also play 
a role in the next subsection and in the discussion of 
the anomalous Hall conductivity and related issues in 
Sec lVLCl 



4. Topological insulators and orbital magnetoelectric response 

There has recently been a blossoming of interest in so- 
called topological insulators, i.e., insulators that cannot 
be adiabatically connected to ordinary insulators with- 
out a gap closure 



Chern-Simons magnetoelectric coupling has been evalu- 
ated from first-principles with the help of WFs for both 



2011]. 

to the case in 



topological and ordinary insulators ( |Coh et al. 

A careful generalization of Eq. ( |94[ ) 
which a finite electric field is present has been carried out 



by Malashevich et al. (2010) in the Wannier representa- 



tion using arguments similar to those in Sees. V.A.I and 
|V.C.2| and used to derive a complete expression for the 
orbital magnetoelectric response, of which the topolog- 



ical Chern-Simons term is only one contribution (Essin 



et al. 2010 Malashevich et al. 2010) 



Hasan and Kane (2010) and Hasan VI. WANNIER INTERPOLATION 



and Moore (2011 1 provide excellent reviews of the back- 



ground, current status of this field, and provide references 
into the literature. 

One can distinguish two kinds of topological insula- 
tors. First, insulators having broken time-reversal (T) 
symmetry (e.g., insulating ferromagnets and ferrimag- 
nets) can be classified by an integer "Chern invariant" 
that is proportional to the Brillouin-zone integral of the 
Berry curvature T n \i summed over occupied bands n. 
Ordinary insulators are characterized by a zero value of 
the invariant. An insulator with a non-zero value would 
behave like an integer quantum Hall system, but without 
the need for an external magnetic field; such systems are 
usually denoted as "quantum anomalous Hall" (QAH) 
insulators. While no examples are known to occur in na- 
ture, tight-binding models exhibiting such a behavior are 



Localized Wannier functions are often introduced in 
textbooks as a formally exact localized basis spanning a 
band, or a group of bands, and their existence provides 
a rigorous justification for the tight-binding (TB) inter- 



1980) 



polation method (Ashcroft and Mermin, 1976 Harrison 



not hard to construct (Haldane 1988). It can be shown 



that a Wannier representation is not possible for a QAH 



In this section we explore the ways in which WFs can 
be used as an exact or very accurate TB basis, allowing 
to perform, very efficiently and accurately, a number of 
operations on top of a conventional first-principles calcu- 
lation. The applications of this "Wannier interpolation" 
technique range from simple band-structure plots to the 
evaluation of various physical quantities as BZ integrals. 
The method is particularly useful in situations where a 
very fine sampling of the BZ is required to converge the 
quantity of interest. This is often the case for metals, as 
the presence of a Fermi surface introduces sharp discon- 
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F(R) 
* ♦ * * * 




FIG. 25 (Color online) Schematic overview of the Wannier interpolation procedure. The left panel shows the BZ mesh q 
used in the first-principles calculation, where the quantity of interest /(q) is explicitly calculated. The Wannier-transformed 
quantity F(R) is strongly localized near the origin of the equivalent supercell, shown in the middle panel with WFs at the 
lattice points. The right panel shows a dense mesh of interpolation points k in the BZ, where the quantity /(k) is evaluated 
at low cost starting from F(JV). 



tinuities in k-space. 

The Wannier interpolation procedure is depicted 
schematically in Fig. [25| The actual first-principles cal- 
culation is carried out on a relatively coarse uniform 
reciprocal-space mesh q (left panel), where the quan- 
tity of interest /(q) is calculated from the Bloch eigen- 
states. The states in the selected bands are then trans- 
formed into WFs, and /(q) is transformed accordingly 
into F (R) in the Wannier representation (middle panel). 
By virtue of the spatial localization of the WFs, -F(R) 
decays rapidly with |R|. Starting from this short-range 
real-space representation, the quantity / can now be ac- 
curately interpolated onto an arbitrary point k in recipro- 
cal space by carrying out an inverse transformation (right 
panel). This procedure will succeed in capturing varia- 
tions in /(k) over reciprocal lengths smaller than the 
first-principles mesh spacing Aq, provided that the lin- 
ear dimensions L = 2n/Aq of the equivalent supercell are 
large compared to the decay length of the WFs. 



A. Band-structure interpolation 

The simplest application of Wannier interpolation is 
to generate band-structure plots. We shall describe the 
procedure in some detail, as the same concepts and no- 
tations will reappear in the more advanced applications 
to follow. 

From the WFs spanning a group of J bands, a set of 
Bloch-like states can be constructed using Eq. Q , which 
we repeat here with a slightly different notation, 



£ elk ' R i 

R 



Rn) 



(n 



1, 



J), (97) 



where the conventions of Eqs. ( 12fl3 ) have been adopted. 
This has the same form as the Bloch-sum formula in 
tight-binding theory, with the WFs playing the role of 
the atomic orbitals. The superscript W serves as a re- 
minder that the states \tp^) are generally not eigenstates 



of the Hamiltonian^] We shall say that they belong to 
the Wannier gauge. 

At a given k, the Hamiltonian matrix elements in the 
space of the J bands is represented in the Wannier gauge 
by the matrix 



W I rr|„ /,W \ _ \ " „ik-R, 



R 



On\H\Bm). (98) 



In general this is a non-diagonal matrix in the band-like 
indices, and the interpolated eigenenergies are obtained 
by diagonalization, 



k.ttm 



(99) 



In the following, it will be useful to view the uni- 
tary matrices as transforming between the Wannier 
gauge on the one hand, and the Hamiltonian (H) gauge 
(in which the projected Hamiltonian is diagonal) on the 
otherj^] From this point forward we adopt a condensed 
notation in which band indices are no longer written ex- 
plicitly, so that, for example, H^ nm = {ip^nWl^km) is 
now written as — (ip^\H and matrix multi- 
plications are implicit. Then Eq. (99) implies that the 



transformation law for the Bloch states is 



(100) 



If we insert into Eqs. (98) and (99) a wavevector be- 



longing to the first-principles grid, we simply recover the 



15 In Ch. Illthe rotated Bloch states IV'^k) were denoted by |i/> n k)) 
see Eq. Hfe). 

16 The unitary matrices C/k are related to, but not the same as, 
the matrices C/' k ' introduced in Eq. |8|. The latter are obtained 
as described in Sees. [ILB] and [H.C| In the present terminology, 
they transform from the Hamiltonian to the Wannier gauge on 
the mesh used in the first-principles calculation. Instead, f/k 
transforms from the Wannier to the Hamiltonian gauge on the 
interpolation mesh. That is, the matrix is essentially an 
interpolation of the matrix [C/' k 'l . 
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first-principles eigenvalues e„k, while for arbitrary k the 
resulting e„k interpolate smoothly between the values on 
the grid. (This is strictly true only for an isolated group 
of bands. When using disentanglement, the interpolated 
bands can deviate from the first-principles ones outside 
the inner energy window, as discussed in Sec. II. I in con- 
nection with Fig. [5]) 

Once the matrices (0|Zf|R) have been tabulated, the 
band structure can be calculated very inexpensively 



by Fourier transforming [Eq. (98 1] and diagonalizing 
[Eq. (99)] matrices of rank J. Note that J, the num- 



ber of WFs per cell, is typically much smaller than the 
number of basis functions (e.g., plane waves) used in the 
first-principles calculation. 

In practice the required matrix elements are obtained 



by inverting Eq. ( 98 ) over the first-principles grid, 



(0|i?|R) = l^e-^ R «|^K) 
q 



(101) 



Here iV is the number of grid points, Eq is the diagonal 
matrix of first-principles eigenenergies, and V q is the ma- 
trix defined in Eq. ( 50 ) , which converts the J7q ah initio 
eigenstates into the J < Jq Wannier-gauge Bloch states, 



(102) 



The strategy outlined above ( Souza et al. 2001 1 is es- 
sentially the Slater-Koster interpolation method. How- 
ever, while the Hamiltonian matrix elements in the local- 
ized basis are treated as adjustable parameters in empir- 
ical TB, they are calculated from first-principles here. 
A similar interpolation strategy is widely used to ob- 
tain phonon dispersions starting from the interatomic 
force constants calculated with density-functional per- 
turbation theory (Baroni et al.. 2001). We shall return 



to this analogy between phonons and tight-binding elec- 



trons (Martin 2004) when describing the interpolation 



of the electron-phonon matrix elements in Sec. |VI.D| 

Wannier band-structure interpolation is extremely ac- 
curate. By virtue of the exponential localization of the 
WFs within the periodic supercell (see Footnote [2| , the 
magnitude of the matrix elements (0|_ff|R) decreases 
rapidly with |R|, and this exponential localization is pre- 
served even in the case of metals, provided a smooth 
subspace has been disentangled. As the number of lat- 



tice vectors included in the summation in Eq. ( 98 ) equals 



the number of first-principles mesh points, beyond a cer- 
tain mesh density the error incurred decreases exponen- 
tially (Yates et al.. 2007). In the following we will il- 



lustrate the method with a few representative examples 
selected from the literature. 




FIG. 26 Band structure of bcc Fe with spin-orbit coupling 
included. Solid lines: original band structure from a con- 
ventional first-principles calculation. Dotted lines: Wannier- 
interpolated band structure. The zero of energy is the Fermi 
level. From|Wang et al\§MX>\. 



1. Spin-orbit-coupled bands of bcc Fe 

As a first application, we consider the relativistic band 
structure of bcc Fe. Because of the spin-orbit interaction, 
the spin density is not perfectly collinear, and the Bloch 
eigenstates are spinors. As mentioned in Sec. |V.C.1[ 
spinor WFs can be constructed via a trivial extension of 
the procedure described in Sec. |TT] for the non-magnetic 
(spinless) case. It is also possible to further modify the 
wannierization procedure so as to produce two separate 
subsets of spinor WFs: one with a predominantly spin-up 
character, and the other with a predominantly spin-down 



character ( Wang et al. 2006 ) 



Using this modified procedure, a set of nine disentan- 
gled WFs per spin channel was obtained for bcc Fe by 
Wang et al. (2006), consisting of three t 2g d-like atom- 
centered WFs and six sp 3 d 2 -\ike hybrids pointing along 
the cubic directions. An inner energy window was cho- 
sen as indicated in Fig. 



26 



so that these 18 WFs describe 
exactly all the the occupied valence states, as well as the 
empty states up to E win , which was set at approximately 
18 eV above the Fermi level. 

The interpolated bands obtained using an 8 x 8 x 8 q- 
grid in the full BZ are shown as dashed lines in Fig. [26] 
The comparison with the first-principles bands (solid 
lines), reveals essentially perfect agreement within the in- 
ner window. This is even more evident in Fig. 27 where 
we zoom in on the interpolated band structure near the 
Fermi level along T-H, and color-code it according to the 
spin-projection along the quantization axis. The vertical 
dotted lines indicate points on the q-mesh. For compar- 
ison, we show as open circles the eigenvalues calculated 
directly from first-principles around a weak spin-orbit- 
induced avoided crossing between two bands of opposite 
spin. It is apparent that the interpolation procedure suc- 
ceeds in resolving details of the true band structure on a 
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FIG. 27 (Color online) Wannier-interpolated relativistic band 
structure of ferromagnetic bcc Fe along F-H. The bands are 
color-coded according to the expectation value of S z : red for 
spin up and blue for spin down. The energies in eV are re- 
ferred to the Fermi level. The vertical dashed lines indicate 
k-points on the mesh used in the first-principles calculation 
for constructing the WFs. From |Yates et al.| ( f2007[ ). 



scale much smaller than the spacing between q-points. 
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FIG. 28 Left panel: band structure of a (5,5) carbon nan- 
otube, calculated by Wannier interpolation (solid lines), or 
from a full diagonalization in a planewave basis set (circles). 
The five vertical dashed lines indicate the five k-points corre- 
sponding to the F point in a 100-atom supercell. The middle 
and right panels show the Wannier-based calculation of the 
quantum conductance and the density of states (see Sec. VII), 
with a perfect match of steps and peaks with respect to the 



exact band structure. From Lee et al. (2005) 



2. Band structure of a metallic carbon nanotube 



As a second example, we consider Wannier interpola- 
tion in large systems (such as nanostructures), that are 
often sampled only at the zone center. We consider here 
a (5,5) carbon nanotube, studied in a 100-atom super- 
cell (i.e. five times the primitive unit cell) and with T- 
point sampling. The system is metallic, and the disen- 
tanglement procedure is used to generate well-localized 
WFs, resulting in either bond-centered combinations of 
sp 2 atomic orbitals, or atom-centered p z orbitals. The 
energy bands at any other k-points are calculated by di- 



agonalizing Eq. ( 98 ) , noting that the size of the supercell 



has been chosen so that the Hamiltonian matrix elements 
on the right-hand-side of this equation are non negligible 
only for WFs up to neighboring supercells on either 
side of R = 0. Fig. [28] shows as solid lines the interpo- 
lated bands, unfolded onto the 20-atom primitive cell. 
Even if with this sampling the system has a pseudogap 
of 2eV, the metallic character of the bands is perfectly 
reproduced, and these are in excellent agreement with 
the bands calculated directly on the primitive cell by di- 
rect diagonalization of the Kohn-Sham Hamiltonian in 
the full plane- wave basis set (open circles). The vertical 
dashed lines indicate the equivalent first-principles mesh 
obtained by unfolding the T-point ^\ 



17 When r sampling is used, special care should be used in calculat- 
ing matrix elements between WFs, since the center of a periodic 
image of, e.g., the ket could be closer to the bra that the actual 
state considered. Similar considerations apply for transport cal- 
culations, and might require calculation of the matrix elements 
in real-space | |Lee| 2006| . 



3. GW quasiparticle bands 

In the two examples above the WFs were generated 
from Kohn-Sham Bloch functions, and the eigenvalues 



used in Eq. (101 1 were the corresponding Kohn-Sham 



eigenvalues. Many of the deficiencies of the Kohn-Sham 
energy bands, such as the underestimation of the energy 
gaps of insulators and semiconductors, can be corrected 
using many-body perturbation theory in the form of the 



GW approximation (for a review, see Aryasetiawan and 
Gunnarsson (1998)). 



One practical difficulty in generating GW band struc- 
ture plots is that the evaluation of the quasiparticle (QP) 
corrections to the eigenenergies along different symmetry 
lines in the BZ is computationally very demanding. At 
variance with the DFT formalism, where the eigenener- 
gies at an arbitrary k can be found starting from the 
self-consistent charge density, the evaluation of the QP 
corrections at a given k requires a knowledge of the 
Kohn-Sham eigenenergies and wavefunctions on a ho- 
mogeneous grid of points containing the wavevector of 
interest. What is often done instead is to perform the 
GW calculation at selected k-points only, and then de- 
duce a "scissors correction," i.e., a constant shift to be 
applied to the conduction-band Kohn-Sham eigenvalues 
elsewhere in the Brillouin zone. 

As already mentioned briefly in Sec 



II. J Hamann 



and Vanderbilt (2009) proposed using Wannier interpo- 



lation to determine the GW QP bands very efficiently 
and accurately at arbitrary points in the BZ. The wan- 
nierization and interpolation procedures are identical to 
the DFT case. The only difference is that the start- 
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FIG. 29 (Color online) Wannier-interpolated upper valence 
and lower conduction bands of SrTiOs from LDA (solid red) 
and GW (dashed blue) calculations. The open circles at the 
symmetry points denote the exact GW results taken directly 
from the first-principles calculation. From [Hamann a nd Van^| 
|derbilt|p009] |. 



ing eigenenergies and overlaps matrices over the uniform 
first-principles mesh are now calculated at the GW level. 
(In the simplest GqWq approximation, where only the 
eigenenergies, not the eigenfunctions, are corrected, the 
wannierization is done at the DFT level, and the resulting 
transformation matrices are then applied to the corrected 
QP eigenenergies.) 

Figure[29|shows a comparison between the interpolated 
GW (dashed lines) and DFT-LDA (solid lines) bands of 
SrTi0 3 ( |Hamann and Vanderbilt] |2009| ) . Note that the 
dashed lines pass through the open circles at the symme- 
try points, which denote exact (non-interpolated) GW 
results. 

Among the recent applications of the GW+Wannicr 
method, we mention the study of the energy bands of zir- 



con and hafnon (Shaltaf et al. 2009), and a detailed com- 



parative study between the DFT-LDA, scissors-shifted, 



and QP GqWq bands of Si and Ge nanowires (Peelaers 



et al. 2011). In the latter study the authors found that 
the simple scissors correction to the DFT-LDA bands 
is accurate near the T point only, and only for bands 
close to the highest valence and lowest conduction band. 



Kioupakis et al. 



(2010) used the method to elucidate 
the mechanisms responsible for free-carrier absorption in 
GaN and in the Ino.25Gao.75N alloy. [Yazyev et aL| ( [2012[ ) 
investigated the quasiparticle effects on the band struc- 
ture of the topolog ical insulators Bi2Se3 and Bi2Te3, and 
Aberg et al. ( 2012| studied in detail the electronic struc- 
ture of LaBr 3 . 



4. Surface bands of topological insulators 

Topological insulators (TIs) were briefly discussed in 



non-magnetic variety, the Z^-odd TIs. The recent flurry 
of activity on this class of materials has been sustained 
in part by the experimental confirmation of the i?2-odd 
character of certain quantum-well structures and of a 
rapidly increasing number of bulk crystals. 

In the case of 3D TIs, the clearest experimental sig- 
nature of the i?2-odd character is at present provided 
by ARPES measurements of the surface electron bands. 
If time-reversal symmetry is preserved at the surface, 
-Z 2 -odd materials possess topologically-protected surface 
states which straddle the bulk gap, crossing the Fermi 
level an odd number of times. These surface states 
are doubly-degenerate at the time-reversal-invariant mo- 
menta of the surface BZ, and in the vicinity thereof they 
disperse linearly, forming Dirac cones. 

First-principles calculations of the surface states for 
known and candidate TI materials are obviously of 
great interest for comparing with ARPES measurements. 
While it is possible to carry out a direct first-principles 
calculation for a thick slab in order to study the topolog- 
ically protected surface states, as done by |Yazyev et al.\ 
(2010), such an approach is computationally expensive. 



Zhang et al. (2010) used a simplified Wannier-based ap- 
proach which succeeds in capturing the essential features 
of the topological surface states at a greatly reduced com- 
putational cost. Their procedure is as follows. First, an 
inexpensive calculation is done for the bulk crystal with- 
out spin-orbit interaction. Next, disentangled WFs span- 
ning the upper valence and low-lying conduction bands 
are generated, and the corresponding TB Hamiltonian 
matrix is constructed. The TB Hamiltonian is then aug- 
mented with spin-orbit couplings AL - S, where A is taken 
from atomic data; this is possible because the WFs have 
been constructed so as to have specified p-like charac- 
ters. The augmented TB parameters are then used to 
construct sufficiently thick free-standing "tight-binding 
slabs" by a simple truncation of the effective TB model, 
and the dispersion relation is efficiently calculated by in- 
terpolation as a function of the wavevector ky in the sur- 
face BZ. 

It should be noted that this approach contains no 
surface-specific information, being based exclusively on 
the bulk WFs. Even if its accuracy is questionable, how- 
ever, this method is useful for illustrating the "topolog- 
ically protected" surface states that arise as a manifes- 



tation of the bulk electronic structure (Hasan and Kane 



2010) 



Sec. V.C.4 (see Hasan and Kane (20101 and Hasan and 
Moore ( 2011[ ) for useful reviews). Here we focus on the 



Instead of applying the naive truncation, it is possible 
to refine the procedure so as to incorporate the changes 
to the TB parameters near the surface. To do so, the 
bulk calculation must now be complemented by a calcu- 
lation on a thin slab, again followed by wannierization. 
Upon aligning the on-site energies in the interior of this 
slab to the bulk values, the changes to the TB param- 
eters near the surface can be inferred. However, [Zhang 
et al. (2011b I found that the topological surface states 
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FIG. 30 Wannier- interpolated energy bands of a free-standing 
(111) slab c ontaining 25 quintuple layers of Bi2Se3 (Zhang 



et al. 



2010) 



plotted along the r-K line in the surface BZ. 
A pair of topologically-protected surface bands can be seen 
emerging from the dense set of projected valence and con- 
duction bulk-like bands and cr ossing at the time- reversal- 
invariant point F. Adapted from Zhang et al. (20101. 
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FIG. 31 (Color online) Surface density-of-states (SDOS) of 
a semi-infinite crystal of Sb2Te3 terminated with a [111] sur- 
face ( Zhang et al.\ |2009a| . Warmer colors represent a higher 
SDOS. The surface states can be clearly seen around Y as red 
lines dispersing in the bulk gap. From Zhang et al. ( 2009a I. 



are essentially the same with and without this surface 
correction. 

The truncated-slab approach was applied by | Zhang 
et al. (2010) to the stoichiometric three-dimensional TIs 



Sb 2 Te 3 , Bi 2 Se 3 , and Bi 2 Te 3 . The calculations on Bi 2 Sc 3 
revealed the existence of a single Dirac cone at the T 
point as shown in Fig. 30 in agreement with ARPES 



measurements (Xia et al. 2009) 



An alternative strategy for calculating the surface 



bands was used earlier by the same authors ( Zhang et al. 



2009a). Instead of explicitly diagonalizing the Wannier- 
based Hamiltonian H(kn) of a thick slab, the Green's 
function for the semi-infinite crystal as a function of 



atomic plane is obtained via iterative methods (Lopez 



Sancho 


et al. 


et al. 


( 


2005 


)• 



1984 1985) 



using the approach of Lee 
Here the localized Wannier representa- 
tion is used to break down the semi-infinite crystal into 
a stack of "principal layers" consisting of a number of 
atomic planes, such that only nearest-neighbor interac- 



tions between principal layers exist (see Ch. VII for more 
details). 

Within each principal layer one forms, starting from 
the fully-localized WFs, a set of hybrid WFs which are 
extended (Bloch-like) along the surface but remain local- 
ized (Wannier-like) in the surface-normal direction (see 



Sees. II. H and V.B.2). This is achieved by carrying out 



a partial Bloch sum over the in-plane lattice vectors, 

l,nk||) = > ; e lK ii-"n|Rn), (103) 



V e lk ii' R n 

Ri 



where I labels the principal layer, kn is the in-plane 
wavevector, and R|| is the in-plane component of R. The 
matrix elements of the Green's function in this basis are 



G 



w (k|.,e) = <k|,/r 



1 



'e-H 



|k||Z'n'). 



(104) 



The nearest-neighbor coupling between principal layers 
means that for each kn the Hamiltonian has a block tri- 
diagonal form (the dependence of the Hamiltonian ma- 
trix on ky is given by the usual Fourier sum expression). 
This feature can be exploited to calculate the diagonal el- 
ements of the Green's function matrix very efficiently us- 
ing iterative schemes (Lopez-Sancho et al. 1984 1985) ^] 
From these, the density of states (DOS) projected onto 



a given atomic plane V can be obtained) Grosso and Par 



ravicim 



2000) as 



iVf (ku,e) = -ilm£ GJJ»(k,|,c + 



(105) 



where the sum over n is restricted to the orbitals ascribed 
to the chosen plane and ry is a positive infinitesimal. 

The projection of the DOS onto the outermost atomic 
plane is shown in Fig. [31] as a function of energy e and 
momentum ky for the (111) surface of Sb 2 Te 3 . The same 
method has been used to find the dispersion of the surface 
bands in the TI alloy Bii-^Sb^ (Zhang et al. 2009b) and 



in ternary compounds with a honeycomb lattice (Zhang 



et al. 2011b). 



B. Band derivatives 

The first and second derivatives of the energy eigenval- 
ues with respect to k (band velocities and inverse effective 
masses) appear in a variety of contexts, such as the cal- 



culation of transport coefficients ( |Ashcroft and Mermin 



18 A pedagogical discussion, where a continued-fractions expansion 
is used to evaluate the Green's function of a semi-infinite linear 
chain with nearest-neighbor interactions, is given by |Grosso and| 

| ParravicTnT| ( |2000} . 
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1976 Grosso and Parravicini 20001. There is therefore 



considerable interest in developing simple and accurate 
procedures for extracting these parameters from a first- 
principles band structure calculation. 

A direct numerical differentiation of the eigenenergies 
calculated on a grid is cumbersome and becomes unre- 
liable near band crossings. It is also very expensive if 
a Brillouin zone integration is to be carried out, as in 
transport calculations. A number of efficient interpola- 
tion schemes, such as the method implemented in the 



Previous first-principles calculations of a xy using vari- 
ous interpolation schemes encountered difficulties for ma- 
terials such as Pd, where band crossings occur at the 
Fermi level (Uehara and Tse, 2000). A Wannier-based 



calculation free from such numerical instabilities was car- 



ried out by Yates et al. (2007), who obtained carefully- 



converged values for a xy in Pd and other cubic metals. 

A more general formalism to calculate the electrical 
conductivity tensor in the presence of a uniform mag- 
netic field involves integrating the equations of motion of 
a wavepacket under the field to find its trajectory on the 
been developed for this purpose, but they are still prone Fermi surface (Ashcroft and Mermin 19761. A numer- 



BoltzTraP package (Madsen and Singh 2006), have 



to numerical instabilities near band degeneracies (Uehara 



and Tse 2000 ) . Such instabilities can be avoided by using 



ical implementation of this approach starting from the 
Wannier-interpolated first-principles bands was carried 



a tight-binding parametrization to fit the first-principles out by Liu et al. ( 2009 1 . This formalism is not restricted 



band structure (Mazin et al. 



As shown by Graf and Vogl 



2000 



Schulz et al. 



1992). to cubic crystals, and the authors used it to calculate the 



(19921 and Boykin (1995), Hall conductivity of hep Mg (Liu et al. 20091 and the 



both the first and the second derivatives are easily com- 
puted within tight-binding methods, even in the presence 
of band degeneracies, and the same can be done in a first- 
principles context using WFs. 

Let us illustrate the procedure by calculating the band 
gradient away from points of degeneracy; the treatment 



of degeneracies and second derivatives is given in Yates 
et al. (2007). The first step is to take analytically the 



derivative d a = d/dk a of Eq. (98) 



H% a =d a H? = Y i e ik - 1l iR a (0\H\IL). (106) 



R 



The actual band gradients are given by the diagonal ele- 
ments of the rotated matrix, 



ft 



(107) 



where is the same unitary matrix as in Eq. ( 99 ) 



It is instructive to view the columns of as or- 
thonormal state vectors in the J-dimensional "tight- 
binding space" defined by the WFs. According to 



Eq. (99) the n-th column vector, which we shall denote 
by ||^nk)), satisfies the eigenvalue equation -f/^H^nk)) = 
£nk||0nk))- Armed with this insight, we now recognize 



in Eq. (107) the Hellmann-Feynman result d a e n ^ — 

«^nk||fl«Wl^k»- 



1. Application to transport coefficients 

Within the semiclassical theory of transport, the elec- 
trical and thermal conductivities of metals and doped 
semiconductors can be calculated from a knowledge of 
the band derivatives and relaxation times r n k on the 
Fermi surface. An example is the low-field Hall conduc- 
tivity a xy of non-magnetic cubic metals, which in the con- 
stant relaxation-time approximation is independent of r 
and takes the form of a Fermi-surface integral containing 



the first and second band derivatives (Hurd 1972) 



magnetoconductivity of MgB 2 ( Yang et al. 2008 ) 



Wannier interpolation has also been used to deter- 
mine the Seebeck coefficient in hole-doped LaRhOa 
and CuRh02 (Usui et al] |2009|), in electron-doped 
SrTi0 3 (lUsui et al. 



et al. 2011 ) 



2010), in SiGe nanowires (Shelley 



and Mostofi 20111, and in ternary skutterudites (Volja 



C. Berry curvature and anomalous Hall conductivity 



The velocity matrix elements between Bloch eigen- 



states take the form (Blount 1962) 



(V'nkl^al^rnk) = S nm d a e nk - i(e mk - e„ k ) [-Ak,a] nTn > 

(108) 

where 



[A 



k.aj 



mk/ 



(109) 



is the matrix generalization of the Berry connection of 



Eq. (95) 



In the examples discussed in the previous section the 
static transport coefficients could be calculated from the 



first term in Eq. ( 108 1 , the intraband velocity. The sec- 



ond term describes vertical intcrband transitions, which 
dominate the optical spectrum of crystals over a wide 
frequency range. Interestingly, under certain conditions, 
virtual interband transitions also contribute to the dc 
Hall conductivity. This so-called anomalous Hall effect 
occurs in ferromagnets from the combination of exchange 
splitting and spin-orbit interaction. For a recent review, 
see 



Nagaosa et al. (20101. 



In the same way that WFs proved helpful for evaluat- 
ing <9 k enk, they can be useful for calculating quantities 
containing k-derivatives of the cell-periodic Bloch states, 



such as the Berry connection of Eq. (109). A number of 



properties are naturally expressed in this form. In addi- 
tion to the interband optical conductivity and the anoma- 
lous Hall conductivity (AHC), other examples include the 



42 



electric polarization (Sec. V. A ) as well as the orbital mag- 
netization and magnetoelectric coupling (Sec. |V.C[). 

Let us focus on the Berry curvature T n ^ [Eq. (|96[)], a 
quantity with profound effects on the dynamics of elec- 



trons in crystals ( |Xiao etal. 20101. T u \l can be nonzero 
if either spatial inversion or time-reversal symmetries are 
broken in the crystal, and when present acts as a kind 
of "magnetic field" in k-space, with the Berry connec- 
tion playing the role of the vector potential. This 
effective field gives rise to a Hall effect in ferromagnets 
even in the absence of an actual applied B-field (hence 
the name "anomalous" ) . The AHC is given bjj^] 



,T AH 



BZ 



dk 



Qtot 



(110) 



where Q k °^ = J2 n /«k^nk,a/3 (/nk is the Fermi-Dirac 
distribution function), and we have rewritten the pseu- 
dovector JF„k as an antisymmetric tensor. 

The interband character of the intrinsic AHC can be 
seen by using k • p perturbation theory to write the k- 



derivatives in Eq. (96), leading to the "sum-over-states" 
formula 



(~)tot 

"a/3 



{fn fm) 



{4>n\ftv a \lj) m )(lf) m \hVp\ fa 



(111) 



The AHC of bec Fe and SrRu03 was evaluated from 



the previous two equations by Yao et al. (2004) and 



Fang et al. (2003) respectively. These pioneering first- 



principles calculations revealed that in ferromagnetic 
metals the Berry curvature displays strong and rapid 
variations in k-space (Fig. 32). As a result, an ultra- 



dense BZ mesh containing millions of k-points is often 
needed in order to converge the calculation. This is the 
kind of situation where the use of Wannier interpolation 
can be most beneficial. 

The strategy for interpolating the Berry curvature is 
similar to that used in Sec. |VI.B| for the band gradient. 
One first evaluates certain objects in the Wannier gauge 
using Bloch sums, and then transform to the Hamilto- 
nian gauge. Because the gauge transformation mixes the 
bands, it is convenient to introduce a generalization of 



Eq. ( 96 ) having two band indices instead of one. To this 



end we start from Eq. ( 109 ) and define the matrices 



Q.a/3 = d a Ap - dpA a = i(d a u\dpu) - i(dpu\d a u), (112) 

where every object in this expression should consistently 
carry either an H or W label. Provided that the cho- 
sen WFs correctly span all occupied states, the inte- 
grand of Eq. (110) can now be expressed as O^j = 
V J f fi H 



19 Equation | |110| gives the so-called "intrinsic" contribution to the 
AHC, while the measured effect also contains "extrinsic" con- 
tributions associated with scattering from impurities l |Nagaosa| 
et aZ.|[20Toj >. 




r h PNr H N r pn 

FIG. 32 Upper panel: band structure near the Fermi level of 
bec Fe with the spontaneous magnetization along z. Lower 
panel: the Berry curvature summed over the occupied bands 
[Eq. | 111 |], plotted along the same symmetry lines. The sharp 
spikes occur when two spin-orbit-coupled bands are separated 
by a small energy across the Fermi level, producing a reso- 
nance enhancement. Adapted from | Yao et a/. | ([2004). 



A useful expression for can be obtained with the 
help of the gauge-transformation law for the Bloch states, 



- i„,w\ 



)C/ k [Eq. (100)]. Differentiating both sides 



(112) 



with respect to k a and then inserting into Eq. 
yields, after a few manipulations, 



n%=n af} - [D a ,A ] + [Dp,A a ] -i[D a ,Dp], (113) 

where D a = Wd a U, and A a , Q, a p are related to the 
connection and curvature matrices in the Wannier gauge 
through the definition G\ = U^O^U-^. Using the band- 



diagonal elements of Eq. (113) in the expression for Q °g 
eventually leads to 



^ ^ ifm fn^)(D anm Ap mn 
rna 



(114) 

This is the desired expression, which in the Wannier in- 
terpolation scheme takes the place of the sum-over-states 



formula. In contrast to Eq. (Ill), note that the summa- 



tions over bands now run over the small set of Wannier- 
projected bands. (Alternatively, it is possible to recast 



Eq. (114) in a manifestly gauge-invariant form such that 
the trace can be carried out directly in the Wannier 



gauge; this formulation was used by Lopez et al. (2012) 



to compute both the AHC and the orbital magnetization 
of ferromagnets.) 



The basic ingredients going into Eq. ( 114 ) are the Wan- 



nier matrix elements of the Hamiltonian and of the posi- 
tion operator. From a knowledge of (0|_ff |R) the energy 
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TABLE I Anomalous Hall conductivity in S/cm of the fer- 
romagnetic transition metals, calculated from first-principles 
with the magnetization along the respective easy axes. The 
first two rows show values obtained using the Wannier in- 
terpolation scheme to either integrate the Berry curvature 
over the Fermi sea (see main text), or to evaluate the Berry 



phases of planar loops around the Fermi surface ( Wang et al. 

20l)7i. Results obtained using the sum-over-states formula 

Eq. Illip, are included for comparison, as well as representa- 
tive experimental values. Adapted from |Wang et al\ ( |2007[ ). 





bec Fe 


fee Ni 


hep Co 


Berry curvature 


753 


-2203 


477 


Berry phase 


750 


-2275 


478 


Sum-over-states 


7530 


-2073^] 


492^] 


Experiment 


1032 


-646 


480 



a |Yao et aE||2004[ |. 

b Y. Yao, private communication. 



gauge-dependent. Nevertheless, the Berry phase ob- 
tained by integrating over a closed loop in k-space, 
ip n = § A„k-dl, is gauge-invariant (Xiao et al. 2010 1. Re- 



calling that T '„k = Vk x A„ k [Eq- ( 96 )J and using Stokes' 
theorem, Eq. (110 1 for the AHC can be recast in terms 



of the Berry phases of Fermi loops on planar slices of 
the Fermi surface. This approach has been implemented 



by Wang et al. (2007), using Wannier interpolation to 
sample efficiently the orbits with the very high density 
required near band-crossings. Table [I] lists values for the 
AHC calculated using both the Berry curvature ( "Fermi 
sea") and Berry-phase ("Fermi surface") approaches. 

It should be possible to devise similar Wannier inter- 
polation strategies for other properties requiring dense 
BZ sampling, such as the magnetic shielding tensors of 
metals ( d' Avezac et al. , 2007 ) . In the following we dis- 



cuss electron-phonon coupling, for which Wannier-based 
methods have already proven to be of great utility. 



eigenvalues and occupation factors, as well as the ma- 
trices U and D a , can b e found using band-structure in- 
terpolation (Sec. VI.Al. The information about A™ and 



f2^|g is instead encoded in the matrix elements (0|r|R) 



as can be seen by inverting Eq. ( 23 ) 



A W = £ ei k.R (0 | ra | R) 



(115) 



R 



As for fijjg, according to Eq. ( 112 ) it is given by the curl 
of this expression, which can be taken analytically. 
The strategy outlined above was demonstrated by 



D. Electron-phonon coupling 



The electron-phonon interaction (Grimvall 1981) 



plays a key role in a number of phenomena, from su- 
perconductivity to the resistivity of metals and the tem- 
perature dependence of the optical spectra of semicon- 
ductors. The matrix element for scattering an electron 
from state ?/> n k to state ^> mj k+q while absorbing a phonon 
qy is proportional to the electron-phonon vertex 



Qv.mn 



(116) 



Wang et al. (| 2006[ ) in calculating the AHC of bec Fe, Here d^V is the derivative of the self-consistent potential 

with respect to the amplitude of the phonon with branch 
index v and momentum q. Evaluating this vertex is a key 
task for a first-principles treatment of electron-phonon 
couplings. 

State-of-the-art calculations using first-principles 
linear-response techniques (Baroni et al. 20011 have 



using the spinor WFs of Sec. |VI.A.1| Both the k-space 
distribution of the Berry curvature and the integrated 
AHC were found to be in excellent agreement with the 



sum-over-states calculation of Yao et al. ( 2004 I 



Table [T] lists the AHC of the ferromagnetic transition 
metal elements, calculated with the magnetization along 
the respective easy axes. The magnetic anisotropy of the 
AHC was investigated by Roman et al. (2009). While 



been successfully applied to a number of problems, 



starting with the works of Savrasov et al. (1994) and 



the AHC of the cubic metals Fe and Ni is fairly isotropic, 
that of hep Co was found to decrease by a factor of four as 
the magnetization is rotated from the c-axis to the basal 
plane. The Wannier method has also been used to calcu- 



Mauri et al. (1996b I, who used respectively the LMTO 
and planewave pseudopotential methods. The cost of 



evaluating Eq. (116) from first-principles over a large 



late the AHC in FePt and FePd ordered alloys ( Seeman 



et al. 2009 Zhang et al. 2011a), and the spin-Hall con- 



ductivity in a number of metals (Freimuth et al. 2010). 



As already mentioned, for certain applications the 
Berry connection matrix [Eq. ( 109 )] is the object of direct 
interest. The interpolation procedure described above 
can be directly applied to the off-diagonal elements de- 
scribing vertical interband transitions, and the magnetic 
circular dichroism spectrum of bec Fe has been deter- 



mined in this way (Yates et al. 2007) 



number of (k, q)-points is quite high, however, and this 
has placed a serious limitation on the scope and accu- 
racy of first-principles techniques for electron-phonon 
problems. 

The similarity between the Wannier interpolation of 
energy bands and the Fourier interpolation of phonon 
dispersions was already noted. It suggests the possibil- 
ity of interpolating the electron-phonon vertex in both 
the electron and the phonon momenta, once Eq. ( |116 ) 
has been calculated on a relatively coarse uniform (k, q)- 
mesh. Different electron-phonon interpolation schemes 



The treatment of the diagonal elements of the Berry 
connection matrix is more subtle, as they are locally 



have been put forth in the literature (Calandra et al. 



2010 Eiguren and Ambrosch-Draxl 2008[ Giustino et al. 
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2007b) In the following we describe the approach first 



developed by Giustino et al. (2007a) and implemented 



in the software package EPW (Noffsinger et al. 2010). 



To begin, let us set the notation for lattice dynam- 



ics (Maradudin and Vosko 19681. We write the instanta- 
neous nuclear positions as R+t s +ur s (£), where R is the 
lattice vector, r s is the equilibrium intracell coordinate 
of ion s = 1, . . . ,p, and UR S (t) denotes the instantaneous 
displacement. The normal modes of vibration take the 
form 



(117) 



The eigenfrequencies aj q „ and mode amplitudes u ql/ 
are obtained by diagonalizing the dynamical matrix 
[l?P h ]™ /3 , where a and (3 denote spatial directions. It is 
expedient to introduce composite indices fj, — (s, a) and 
v = (t,/3), and write for the dynamical matrix. 

With this notation, the eigenvalue equation becomes 



[«e q ] 



(118) 
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FIG. 33 (Color online) Real-space representation of 
(O e |9R pM V|R e ), the electron-phonon vertex in the Wannicr 
basis. The black squares denote the crystal lattice, the blue 
lines the electron Wannier functions |0 e ) and |R) e , and the 
red line the phonon perturbation in the lattice Wannier rep- 
resentation, dn. p V(r). Whenever two of these functions are 
centered on distant unit cells, the vertex becomes vanishingly 
small. From |Giustino et aT| ( |2007a| >. 



where e q is a 3p x 3p unitary matrix. In analogy with 
the tight-binding eigenvectors 
view the columns of e r 



lk )) of Sec. |VI.B 
-q./ii/ as orthonormal phonon eigen- 
vectors e qiy . They are related to the complex phonon 
amplitudes by u qi/ = (mo/m s ) 1 / 2 e qi/ (mo is a reference 
mass), which we write in matrix form as U^\ L 



we can The object (0 c |9r ^VjRg), the electron-phonon vertex 



Returning to the electron-phonon vertex, Eq. ( 116 1, we 



can now write explicitly the quantity dq U V therein as 



9 q ,F(r) = | ? y(r;{R 



?7 U R S }) 



E e^ R d R ,V(r)U^, 

R,fj, 



(119) 



where ^^(r) is the derivative of the self-consistent po- 
tential with respect to mr s ,q- As will be discussed in 
Sec. |VIII| it is possible to view these single-atom dis- 
placements as maximally-localized "lattice Wannier func- 
tions." With this interpretation in mind we define the 
Wannier-gauge counterpart of <9 qt/ V(r) as 



R 



e iq -%/(r) 



(120) 



related to the "eigenmode-gauge" quantity 9 qi/ V(r) by 
d^V(v) = E d™V(r)U^ u . (121) 



in the Wannier representation, is depicted schematically 
in Fig. 33 Its localization in real-space ensures that 
Eq. (122) smoothly interpolates g w in the electron and 



phonon momenta. Finally, we transform the interpolated 
vertex back to the Hamiltonian/cigenmode gauge, 

^(k',q') = ^w|S q v^|^) 



k'+q 



-[Es^^C]^ (123) 



where we used Eqs. (100) and (121). 

Once the matrix elements (0 e |9R p;i V*|R e ) are known, 
the electron-phonon vertex can be evaluated at arbitrary 
(k', q') from the previous two equations. The procedure 
for evaluating those matrix elements is similar to that 
leading up to Eq. (floi! for (0|iT|R): invert Eq. Jl22l 



over the first-principles mesh, and then use Eqs. ( 102 ) 



and (121) 



The above interpolation scheme has been applied to 
a number of problems, including the estimation of T c 



in superconductors (Giustino et al. 2007b Noffsinger 



et al. 2008 2009); the phonon renormalization of en- 



ergy bands near the Fermi level in graphene ( Park et al. 



2007 ) and copper oxide superconductors ( |Giustino et al. 



2008); the phonon renormalization of the band gap of 
diamond (Giustino et al. 2010); the vibrational life- 



Next we introduce the Wannier-gauge vertex (k, q) = times ( Park et al. 2008 ) and electron linewidths ( Park 
(^k^ql^q^^l^k')' which can be readily interpolated onto 
an arbitrary point (k',q') using Eqs. ([97]) and (120), 



et al. 2009 ) arising from electron-phonon interactions in 



graphene; and the phonon-assisted optical absorption in 



q W 



(k',q')= £ e l ( k '- R ^ q '- R ^(O e |G>R p ^|R e 



(122) 



silicon (Noffsinger et al. 2012) (in this last application 
both the velocity and the electron-phonon matrix ele- 
ments were treated by Wannier interpolation). 



We mention in closing the work of |Calandra et al. 
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(2010), where the linear-response calculation of the de- 
formation potential dq„V is carried out taking into ac- 
count nonadiabatic effects (that is, going beyond the 
usual approximation of static ionic displacements). Us- 
ing the electron-phonon interpolation scheme of|Giustino| 



across length-scales. In particular, Hierse and Stechel 



et al. (2007a I to perform the BZ integrations, and a Wan- 



nier interpolation of the dynamical matrix to capture 
Kohn anomalies, the authors found significant nonadia- 
batic corrections to the phonon frequencies and electron- 
phonon coupling matrix elements in MgB 2 and CaCg- 



VII. WANNIER FUNCTIONS AS BASIS FUNCTIONS 



In Ch. VI we described the use of Wannicr functions 



as a compact tight-binding basis that represents a given 
set of energy bands exactly, and that can be used to cal- 
culate a variety of properties efficiently, and with very 
high accuracy. This is possible because (a) WFs and 
bands are related by unitary transformations, and (b) 
WFs are sufficiently localized that any resulting tight- 
binding representation may be truncated with little loss 
of accuracy. In this Chapter, we review two further gen- 
eral approaches to the use of WFs as optimal and com- 
pact basis functions for electronic-structure calculations 
that exploit their localization and transferability. 

The first category includes methods in which WFs are 
used to go up in the length scale of the simulations, using 
the results of electronic-structure calculations on small 
systems in order to construct accurate models of larger, 
often meso-scale, systems. Examples include using WFs 
to construct tight-binding Hamiltonians for large, struc- 
turally complex nanostructures (in particular for study- 
ing quantum transport properties), to parametrize semi- 
empirical force-fields, and to improve the system-size 
scaling of quantum Monte Carlo (QMC) and GW calcu- 
lations, and the evaluation of exact-exchange integrals. 

The second category includes methods in which WFs 
are used to identify and focus on a desired, physically rel- 
evant subspace of the electronic degrees of freedom that 
is singled out ( "downfolded" ) for further detailed anal- 
ysis or special treatment with a more accurate level of 
electronic structure theory, an approach that is particu- 
larly suited to the study of strongly-correlated electron 
systems, such as materials containing transition metal, 
lanthanoid, or even actinoid ions. 



A. WFs as a basis for large-scale calculations 



Some of the first works on linear-scaling electronic 



structure algorithms (Baroni and Giannozzi 



1992 



and Parrineho] [19921 |Hierse and Stechel[ |1994 



Galli 



Yang 



1991 ) highlighted the connection between locality in elec- 



tronic structure, which underpins linear-scaling algo- 



(1994) considered explicitly two large systems A and B, 



different globally but which have a certain similar local 
feature, such as a particular chemical functionalization 
and its associated local environment, which we will call 
C. They argued that the local electronic structure infor- 
mation associated with C should be similar whether cal- 
culated from system A or system B and, therefore, that 
it should be possible to transfer this information from a 
calculation on A in order to construct a very good ap- 
proximation to the electronic structure in the locality of 
feature C in system B. In this way, large computational 
savings could be made on the self-consistency cycle, en- 
abling larger-scale calculations. 

The units of electronic structure information that 



Hierse and Stechel (1994) used were localized non- 



orthogonal orbitals, optimized in order to variationally 
minimize the total energy of the system. These orbitals 
are also referred to in the literature as "support func- 
tions" (Hernandez and Gillan 1995) or "non-orthogonal 



generalized Wannier functions" (Skylaris et al. 2002). 



1. MLWFs as electronic-structure building blocks 

Since MLWFs encode chemically accurate, local (and 
thus transferable) information, they can act as building 
blocks to construct the electronic structure of very large- 
scale systems (ILee et al] 120051). In this approach the 



Hamiltonian matrix of a large nanostructure for which a 
full, conventional DFT calculation would be intractable, 
is built using first-principles calculations performed on 
smaller, typically periodic fragments. The matrix ele- 
ments in the basis of MLWFs that are obtained from the 
calculations on the fragments can be used to construct 
the entire Hamiltonian matrix of the desired system, with 
the size of the fragments a systematically controllable de- 



terminant of the accuracy of the method (Shelley et al. 



2011). Such an approach has been applied to complex 



systems containing tens of thousands of atoms (Cantele 



et al. 20091 Lee and Marzari, 2006 



Li and Marzaril 2011 



Li et al.\ [20TT] |Shelley and Moston] [20TT] |Shelley et al. 
20111). 



rithms (Sec. III.D ), and the transferability of information 



An issue that arises when combining matrix elements 
from more than one DFT calculation into a single tight- 
binding Hamiltonian is that a common reference poten- 
tial must be defined. For example, consider combining 
the results of a calculation on a perfect bulk material 
and one on the same material with an isolated structural 
defect. In the latter case, the diagonal (on-site) matrix 
elements (w n \H\w n ) in the system with the defect should 
converge to the same value as those in the pristine sys- 
tem as one goes further away from the location of the de- 
fect. With periodic boundary conditions, however, this 
is not guaranteed: the difference between the matrix ele- 
ments in the respective calculations will, in general, tend 
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FIG. 34 Top: Illustration of the lead-conductor-lead geome- 
try used in quantum transport calculations. The conductor 
(C), left lead (L) and right lead (R) are described by Hamil- 
tonians He, Hl and Hr, respectively. The coupling between 
adjacent regions is described by matrices hhc and hen- Bot- 
tom: The leads are split into principal layers corresponding 
to the sub-matrices in Eq. \12&\ . From |Shelley et fflZ.|f2011| ). 



to a non-zero constant, this being the difference in ref- 
erence potential between the two calculations. In other 
words, the difference in value between the on-site matrix 
elements from the bulk calculation and those from the 
calculation with the defect, but far away from the loca- 
tion of the defect, give a direct measure of the potential 
offset between the two calculations. This offset can then 
be used to align all of the diagonal matrix elements from 
one calculation with those of the other, prior to them be- 
ing combined in a tight-binding Hamiltonian for a larger 
nanostructure. The potential alignment approach de- 
scribed above has been used for transport calculations 
(see next subsection), and for the determination of band 
offsets and Schottky barriers ( Singh-Miller 2009 1 , for cal- 



culating formation energies of point defects ( Corsetti and 



Mostofi 2011), and for developing tight-binding models 



of the surfaces of topological insulators (Zhang et al. 



2010) (see Sec. VI.A.4) 



Last, it should be noted that charge self-consistency 
could play an important role when trying to build the 
electronic structure of large or complex nanostructurcs, 
and one might have to resort to more sophisticated ap- 
proaches. As a suggestion, electrostatic consistency could 
be attained solving the Poisson equation for the entire 



system (Leonard and Tersoff 1999), then using the up- 



dated electrostatic potential to shift appropriately the 
diagonal elements of the Hamiltonian. In a more general 
fashion, the electronic charge density and the Hartree and 
exchange-correlation potentials could be updated and op- 
timized self-consistently in a basis of disentangled, frozen 



MLWFs, using a generalized occupation matrix (Marzari 



et al. 1997). 



2. Quantum transport 

A local representation of electronic structure is partic- 
ularly suited to the study of quantum transport, as we 
illustrate here in the case of quasi-one-dimensional sys- 
tems. We consider a system composed of a conductor C 
connected to two semi-infinite leads, R and L, as shown 



in Fig. 34 The conductance Q through a region of in- 
teracting electrons is related to the scattering properties 



of the region itself via the Landauer formula (Landauer 



1970) 



2e 5 



T(E), 



(124) 



where the transmission function T(E) is the probabil- 
ity that an electron with energy E injected at one end 
of the conductor will transmit to the other end. This 
transmission function can be expressed in terms of the 
Green's functions of the conductors and the coupling of 



1981). 



the conductor to the leads (Datta 1995 Fisher and Lee 



T(E) = Tv{T L G r c T R G a c ), 



(125) 



where G^' a ^ are the retarded (r) and advanced (a) 
Green's functions of the conductor, and Ti^ m are func- 
tions that describe the coupling of the conductor to the 
left (L) and right (R) leads. Since G a = {G r )\ we con- 
sider G r only and drop the superscript. 

Expressing the Hamiltonian H of the system in terms 
of a localized, real-space basis set enables it to be par- 
titioned without ambiguity into sub-matrices that cor- 
respond to the individual subsystems. A concept that 



is particularly useful is that of a principal layer (Lee 



and Joannopoulos 



1981) (PL), which is a section of 
lead that is sufficiently long such that (xf \H\x™) — if 
|m — n\ > 2, where H is the Hamiltonian operator of the 
entire system and |x") is the i basis function in the n th 
PL. Truncating the matrix elements of the Hamiltonian 
in this way incurs a small error which is systematically 
controlled by increasing the size of the PL. The Hamil- 
tonian matrix in this basis then takes the block diagonal 



form (see also Fig. 34) 
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, (126) 



where He represents the Hamiltonian matrix of the con- 
ductor region, Hf and Hf are those of each PL of the 
left and right leads, respectively, Hf and H R X are cou- 
plings between adjacent PLs of lead, and h^c and hcR 
give the coupling between the conductor and the leads. 

In order to compute the Green's function of the con- 
ductor one starts from the equation satisfied by the 
Green's function G of the whole system, 



(e - H)G 



(127) 
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where I is the identity matrix, and e 
is an arbitrarily small, real constant. 



E + if], where T) Mostofi 2011), as well as more methodological work on 



improving the description of electron correlations within 



From Eqs. (127) and (126), it can be shown that the 



Green's function of the conductor is then given by (Datta 12007 Ferretti et al. 



1995) 



this formalism (Bonferroni et al. 

" 2005b[ 



2008 Calzolari et al. 



Gc — ( e — He — Sl — Tiii) 



(128) 



where we define S L = h\ c g L h LC and T, R = h RC gnh} RC , 
the self-energy terms due to the semi-infinite leads, and 
9{l,r} — (e — HiL^m) , the surface Green's functions 
of the leads. 

The self-energy terms can be viewed as effective Hamil- 
tonians that arise from the coupling of the conductor 
with the leads. Once the Green's functions of the leads 
are known, the coupling functions r^^j can be easily 
obtained as (Datta 1995) 



The formulation described above relies on a localized 
description of the electronic-structure problem, and it 
should be noted that several approaches to calculating 
electronic transport properties have been developed us- 
ing localized basis sets rather than MLWFs, ranging from 
Gaussians (Hod et al. 2006) to numerical atomic or- 



bitals (Brandbyge et al. 



Rocha et al. 2008 ) 



2002 Markussen et al. 2006 



r 



{L,R} 



{L,R} 



^{L,R}\ 



(129) 



It 



where E^ fl} = (E^Jt 

As for the surface Green's functions g{L,R} of the semi- 
infinite leads, it is well-known that any solid (or surface) 
can be viewed as an infinite (semi-infinite in the case of 
surfaces) stack of principal layers with nearest-neighbor 
interactions (Lee and Joannopoulos 1981). This corre- 



In addition, in the Keldysh formalism one can add 
more complex interaction terms to the self-energies, 
such as electron-phonon or (for molecular conductors) 
electron- vibration interactions (Frederiksen et al. 2007[ ). 
These latter can also be conveniently expressed in a ML- 
WFs representation, and a natural extension of the previ- 
ous quantum-transport formalism to the case of inelastic 
scattering channels due to electron- vibration interactions 



has recently been developed in a MLWFs basis (Kim and 



Marzari 2012). 



sponds to transforming the original system into a lin- 
ear chain of PLs. Within this approach, the matrix el- 



ements of Eq. (127) between layer orbitals will yield a 
set of coupled equations for the Green's functions which 
can be solved using an efficient iterative scheme due to 



Lopez-Sancho et al. (1984 19851. Knowledge of the fi- 



nite Hamiltonian sub-matrices in Eq. (126), therefore, is 



sufficient to calculate the conductance of the open lead- 



conductor- lead system given by Eq. (124). 

There are a number of possibilities for the choice of 
localized basis |x). Early work used model tight-binding 



Hamiltonians (Anantram and Govindan 



1998 Chico 



3. Semi-empirical potentials 

First-principles molecular dynamics simulations of 
large-scale (thousands of atoms) systems for long 
(nanoseconds) time-scales are computationally costly, if 
not intractable. Molecular dynamics simulations with 
empirical interatomic potentials are a feasible alternative 
and there is an ongoing effort in developing potentials 
that are more accurate, more transferable and, therefore, 
more predictive. One approach in this direction is to fit 
the parameters that appear within empirical potentials 
so that they reproduce target properties, such as forces 
and stresses, obtained from accurate DFT calculations 



et al. 1996 Nardelli 1999 Saito et al. 1996 ), but the in- on a large number of atomic configurations (Ercolessi 



creasing sophistication of computational methods meant 
that more realistic first-principles approaches could be 



adopted ( 


Buongiorno Nardelli et al. 


2001 


Fattebert and 


Buongiorno Nardelli 


2003 


). Maximally-' 


ocalized Wan- 



et al. (2004), who studied Al and C chains and a (5,0) 



carbon nanotube with a single Si substitutional defect, 



and by Lee et al. ( 2005 ) , who studied covalent function- 
alizations of metallic nanotubes - capabilities now en- 



coded in the open-source packages Wannier90 (Mostofi 



et al.. 2008) and WanT (Ferretti et al. 2005a). This was 



quickly followed by a number of applications to ever more 
realistic systems, studying transport through molecular 



sen 



junctions (Strange et al. 2008 Thygesen and Jacob 



and Adams 1994 1. In the particular class of ionic con- 
densed matter systems, e.g., first and second row metal 
oxides, it is well-known that the electronic properties of 
an ion can be significantly affected by its coordination 
environment and, therefore, that it is also important to 
include an accurate description of polarization effects in 
any interatomic potential. While simple potentials may 
attempt to account for these many-body effects in an av- 
erage manner, the result is that the potential loses trans- 
ferability and, hence, predictive power. As a result, there 
has been some effort in developing potentials such that 
they are also fitted to information from DFT calculations 
regarding the electron distribution around each ion, in 
particular, dipoles and quadrupoles. 



2005|), decorated carbon nanotubes and nanorib 

" Li 



bons (Cantele et al. 2009 Lee and Marzari 2006 



et al. 2011 Rasuli et al. 2010 ), organic monolayers (Bon- 



ferroni et al. 2008), and silicon nanowires (Shelley and 



In this vein, Aguado et al. (2003a) introduced the use 



of MLWFs in order to calculate dipole and quadrupole 
moments in MgO and used these to construct an in- 
teratomic potential based on the aspherical ion method 
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(AIM) potential ( [Rowley et aTj | 1998 |). In an ionic system 



such as this, MLWFs are found to be well-localized close 
to the ions and each can therefore be associated unam- 
biguously with a particular ion. The dipole moment fi 1 
of each ion / is then calculated as 



(130) 



nel 



where a is a Cartesian component, Zi is the charge of 
ion I at position R/ , f na is the center of the n th MLWF 



and is given by Eq. ( 28 ) or Eq. ( 43 1 , and the factor of 



two accounts for spin degeneracy. 

For the quadrupole moments B a a, the fact that the 
MLWFs are localized within the simulation cell is ex- 
ploited in order to explicitly compute the real-space in- 
tegral 



9 7 

'a/3 



2^ / dv K(r)| 2 [3r^-(r')%3]/2, 



where r 1 = r — R 7 , and the integral is performed over a 
sphere Vj , centered on R 1 , such that the integral of the 
electron density associated with the MLWF within this 
sphere is greater than some threshold. The potential ob- 
tained exhibits good transferability and the method has 
been used to parametrize similar potentials for other al- 



kaline earth oxides (Aguado et al. 2003b), A1 2 3 (Jahn 



etal 2006), and the CaO-MgO-Al 2 3 -Si02 (CMAS) 



mineral system (Jahn and Madden 2007). 



The use of MLWFs for attempting to obtain better in- 
teratomic potentials has not been limited to the solid 
state. In biomolecular simulations, an important fac- 
tor in developing accurate force-fields is having an ac- 
curate description of the electrostatic potential. Starting 
from DFT calculations on isolated molecules, |Sagui et al. 
( 2004 ) partition the electronic charge density into contri- 



butions from individual MLWFs and calculate the multi- 
poles of each using an order-by-order expansion of gauge- 



invariant cumulants ( Resta 1998 
reader is referred to Sagui et al. 



Souza et al. 2000) (the 



(2004) for full details). 



Using fast particle mesh Ewald and multigrid methods, 
these multipoles can then be used to generate the elec- 



trostatic potential. Sagui et al. ( 2004 ) show that higher 



order multipoles, e.g., up to hexadecapole, may be in- 
corporated without computational or numerical difficulty 
and that agreement with the 'exact' potential obtained 
from DFT is very good. The idea of partitioning the 
charge density according to individual MLWFs was also 
employed by |Kirchner and Hutter| ( |2004[ ) in order to 
determine atomic charges in dimethyl sulfoxide, show- 
ing that there can be significant deviations between the 
gaseous and aqueous forms and, therefore, underlining 
the importance of using polarizable force-fields for de- 
scribing solvated systems. 



Finally, we note that recently Rotenberg et al. (2010) 



entirely on a partitioning of the electronic density in 
terms of MLWFs. Their method has been applied suc- 
cessfully to water and molten salts. It remains to be seen 
whether the approach is extensible to more complex or 
anisotropic systems. 



4. Improving system-size scaling 

The localized nature of MLWFs in real-space makes 
them a natural and appealing choice of basis for elec- 
tronic structure calculations as the sparsity exhibited by 
operators that are represented in a localized basis may 
be exploited in order to achieve computational efficien- 
cies and improved scaling of numerical algorithms with 
respect to system size. Recently, therefore, MLWFs have 
been used for this very purpose in a number of contexts 
and we mention, in brief, some of them here. 



In quantum Monte Carlo (QMC) calculations (Foulkes 



(131) et al. 20011, a significant computational effort is ex- 



pended in evaluating the Slater determinant for a given 
electronic configuration of N electrons. These Slater 
determinants are usually constructed from a set of ex- 
tended single-particle states, obtained from, e.g., a DFT 
or Hartree-Fock calculation, represented in a basis set of, 
e.g., extended plane- waves. This gives rise to 0(N 3 ) scal- 
ing of the computational cost of evaluating such a deter- 
minant. Williamson et al. (20011 suggested, instead, to 



use MLWFs that were smoothly truncated to zero beyond 
a certain cut-off radius that is independent of system size. 
This ensures that each electron falls only within the local- 
ization region of a fixed number of MLWFs, thus reduc- 
ing the asymptotic scaling by one factor of N. Further- 
more, by representing the MLWFs in a basis of localized 
spline functions, rather than plane-waves or even Gaus- 
sian functions, the evaluation of each orbital is rendered 
independent of system size, thereby reducing the overall 
cost of computing the determinant of the given config- 
uration to O(N). More recently, rather than truncated 
MLWFs, the use of non-orthogonal orbitals obtained by 



projection ( |Reboredo and Williamson| |2005[ ) or other lo- 
calization criteria (Alfe and Gillan 2004) has also been 
suggested. 



In another development, Wu et al. (20091 use ML- 
WFs in order to compute efficiently Hartree-Fock 
exact-exchange integrals in extended systems. Hybrid 



exchange-and-correlation functionals (Becke 1993) for 



have proposed force-fields whose parametrization is based 



DFT calculations, in which some proportion of Hartree- 
Fock exchange is included in order to alleviate the well- 
known problem of self-interaction that exists in local and 
semi-local functionals such as the local-density approxi- 
mation and its generalized gradient-dependent variants, 
have been used relatively little in extended systems. This 
is predominantly due to the computational cost associ- 
ated with evaluating the exchange integrals between ex- 
tended eigenstates that are represented in a plane-wave 
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basis set. Wu et al. (20091 show that by performing a 



unitary transformation of the eigenstates to a basis of 
MLWFs, and working in real-space in order to exploit 
the fact that spatially distant MLWFs have exponen- 
tially vanishing overlap, the number of such overlaps that 
needs to be calculated scales linearly, in the limit of large 
system-size, with the number of orbitals (as opposed to 
quadratically) , which is a sufficient improvement to en- 
able Car-Parrinello molecular dynamics simulations with 
hybrid functionals. 

Similar ideas that exploit the locality of MLWFs have 
been applied to many-body perturbation theory ap- 
proaches for going beyond DFT and Hartree-Fock cal- 
culations, for example, in the contexts of the GW ap- 



proximation (Umari et at, 2009), the evaluation of local 



correlation in extended systems (Buth et al. 2005 Pisani 



et al. 


2005 


et al. 


2010) 



2005), and the Bethe-Salpeter equation (Sasioglu 
et al. 2010). The improved scaling and efficiency of these 
approaches open the possibility of such calculations on 
larger systems than previously accessible. 

Finally, we note that MLWFs have been used recently 
to compute van der Waals (vdW) interactions in an ap- 



proximate but efficient manner ( Andrinopoulos et al. 



2011 Silvestrelli, 2008 



2009b[). The method i s based 



on an expression due to Andersson et al. ( 1996 1 for the 
vdW energy in terms of pairwise interactions between 
fragments of charge density. MLWFs provide a localized 
decomposition of the electronic charge density of a sys- 
tem and can be used as the basis for computing the vdW 
contribution to the total energy in a post-processing (i.e., 
non-self-consistent) fashion. In order to render tractable 
the necessary multi-dimensional integrals, the true ML- 
WFs of the system are substituted by analytic hydrogenic 
orbitals that have the same centers and spreads as the 
true MLWFs. The approach has been applied to a vari- 
ety of systems, including molecular dimers, molecules ph- 
ysisorbed on surfaces, water clusters and weakly-bound 



solids (Andrinopoulos et al. 


2011 


Espejo et al. 


2012 


Silvestrelli 


2008 2009a b 


Silvestrelli et al. 


2009 


1. Re- 


cently, 


Ambrosetti and Silvestrelli ( 


2012 


) have sug 


gested 



an alternative, simpler formulation that is based on Lon- 
don's expression for the van der Waals energy of two 



interacting atoms (Eisenschitz and London 1930). 



B. WFs as a basis for strongly-correlated systems 

For many strongly-correlated electron problems, the 
essential physics of the system can be explained by con- 
sidering only a subset of the electronic states. A re- 
cent example is understanding the behavior of unconven- 
tional (high-T c ) superconductors, in which a great deal 
of insight can be gained by considering only the ML- 
WFs of p and d character on Cu and O, respectively, 



for cuprates (Sakakibara et al. 2010), and those on As 



and Fe, respectively, for the iron pnictides (Cao et al. 



2008} |Hu andTLu1|20T0{|Kuroki et at] |2008{ [Suzuki et al. 



2011 ). Other strongly-correlated materials for which ML- 



WFs have been used to construct minimal models to 



help understand the physics include manganites ( Kovacik 
and Ederer] 2010), topological insulators (Zhang et al. 



2009a|b ) (see also Sec. VI.A.4), and polyphenylene viny 



lene (PPV), in particular relating to electron- hole exci- 
tations ( |Karabunarliev and Bittner 2003a|b ). 

Below we outline some of the general principles behind 
the construction and solution of such minimal models. 



1. First-principles model Hamiltonians 

Consider a strongly-correlated electron system de- 
scribed by a Hamiltonian of the form 



H = Hn 



int j 



(132) 



where Hq contains the one-particle terms and Hi nt the 
interaction (e.g., Coulomb repulsion) terms. In second- 
quantized notation and expressed in terms of a complete 
tight-binding basis, this may be expressed as 

ij R,R' 

+ E ^' R " R "'sk4R' a ^"^'"' 



ijkl RR'R"R"' 



(133) 



where R labels the correlated site, lower-case indices 
(such as i and j) refer to the orbital and spin degrees 
of freedom, c\ R (cm) creates (annihilates) an electron in 
the orbital wm(r), and h and U contain the matrix el- 
ements of the single-particle and (screened) interaction 
parts of the Hamiltonian, respectively. 

Usually, a complete tight-binding representation of all 
the states of the system is not required, and the essential 
physics can be described by a smaller set of physically 
relevant orbitals. For example, those states close to the 
Fermi level, or those of a particularly symmetry, or those 
localized on specific sites, may be sufficient. In this way, 
the size of the basis used to represent the many-body 



Hamiltonian is greatly reduced, rendering Eq. ( 133 ) more 



amenable to solution (see, e.g., Solovyev (2008)). 

In this spirit, MLWFs obtained from DFT calculations 
have been used as the orbital basis for these minimal 
models^] Advantages of using MLWFs include the fact 
that they can be chosen to span precisely the required 
energy range (using the disentanglement procedure out- 
lined in Sec. II.I.2), and that they naturally incorporate 



20 An alternative approach is to obtain the orbitals via the down- 
folding method, discussed in Ch. Ill II 
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hybridization and bonding appropriate to their local en- 
vironment. 

The single-particle hopping parameters of the model 
Hamiltonian are obtained easily from the matrix ele- 
ments of the DFT Hamiltonian represented in the ba- 
sis of MLWFs, using Eq. (101). The interaction pa- 



rameters of the model Hamiltonian can be calculated, 



for example, from either constrained DFT (Anisimov 


et al.\ '. 


L991; Dcderichs et aZ.||1984; McMahan et al. 1990 


Nakamura et al.\ 20061, within the random phase ap- 


proximation ( Aryasetiawan et al. 2004 Miyake et al. 


2009 5 


Solovyev and Imada 2005| |Springer and Aryase-| 


tiawan 


1998), or by direct calculation of the matrix 



elements of a suitable screened Coulomb interaction be- 



tween, for example, MLWFs (Miyake et al. 2006 Naka- 



mura et al. 2006). It is interesting to note that numer- 
ical evidence suggests that on-site Coulomb interactions 
(both screened and bare), are maximized when calcu- 



lated with a basis of MLWFs ( Miyake and Aryasetiawan 



2008 ) and, therefore, that MLWFs may be an optimally 
localized basis for this purpose. This is perhaps not sur- 
prising given the broad similarities between MLWFs and 
WFs obtained via the Edmiston-Ruedenberg localization 



scheme (Edmiston and Ruedenberg 19631, discussed in 



Sec. |III.A which maximizes the electronic Coulomb self- 
interaction of each orbital. 

Once the parameters of the model have been deter- 
mined, the model Hamiltonian is then said to be "from 
first-principles" , in the sense that the parameters of the 
model are determined from DFT rather than by fitting to 
experiments. The many-body Hamiltonian in the mini- 
mal basis of MLWFs may be then solved by direct di- 
agonalization, or one of a number of other approaches 
that are too numerous to review here but which include, 
for example, the dynamical mean-field theory (DMFT) 
approach. DMFT maps the many-body problem on to 



an Anderson impurity model (Anderson 1961) in which 
on-site correlation is treated non-pcrturbatively and the 
correlated sites are coupled to a self-consistent energy 
bath that represents the rest of the system. The impu- 
rity sites, also known as the "correlated subspaces", are 
defined by localized projector functions and MLWFs are 



a common choice (Lechermann et al. 2006 Trimarchi 



et al. 2008 Weber et al. 2010). In particular, one would 



typically choose orbitals of d or / character associated 
with transition metal, lanthanoid or actinoid ions. The 
Green's function for the impurity site is calculated self- 
consistently, for example, by a numerical functional in- 
tegration (which constitutes the bulk of the computa- 
tion). Further self-consistency with the DFT ground- 
state may also be attained by using the solution to the 
impurity problem to update the electronic density that is 
then used to construct an updated Kohn-Sham potential, 
which determines a new set of eigenstates, MLWFs and, 
hence, model Hamiltonian parameters that can then be 
fed back in to the DMFT cycle. The reader is referred 



to Kotliar et al. (2006) and Held (2007) for further de- 



tails; other examples that use localized Wannier functions 
or generalized tight-binding models to address correlated 



electrons problems can be found in 


Amadon et al. 


(2008); 


Anisimov et al. 


(2005 


); 


Held et al. 


(2008 


); 


Korotin et al. 


(2008 


); Korshunov et al. 


(2005 


); and 


Ku et al. 


(2002 


1. 



2. Self-interaction and DFT + Hubbard U 

In the approaches just described, the results of a DFT 
calculation are used to parametrize the model Hamilto- 
nian of a strongly correlated electron system. In contrast, 



in a DFT+C7 formulation (Anisimov et al. 1993 1991) 



the energy functional is explicitly augmented with a Hub- 
bard term U ( |Hubbard[ |1963[ ) aimed at improving the de- 
scription of strong interactions, such as those associated 
with localized d and / electronic states, and at repair- 
ing the mean-field underestimation of on-site Coulomb 
repulsions. 

In DFT+[7 the Hubbard manifold is defined by a set of 
projectors that are typically atomic-like orbitals of d or / 
character. Localization of this manifold plays a key role, 
since DFT+f/ effectively corrects for the lack of piecewise 



linearity in approximate energy functionals (Cococcioni 



and de Gironcoli 


2005 


Perdew et al. 


1982 


) and thus 


greatly reduces self-interaction errors ( 


Kulik et al. 


2006 


Mori-Sanchez et al. 


2008 


Perdew and Zunger 


19811. 



Since strongly localized orbitals are those that suffer most 
from self-interaction, MLWFs can become an appealing 
choice to define Hubbard manifolds adapted to the chem- 
istry of the local environment. In fact, MLWFs have 



been successfully used as Hubbard projectors (Anisimov 
et al.\ |2007t |Fabris et~aT\ |2005| |Miyake and AryasetT 



awan 



2008), and it has been argued that their shape can 



constitute an additional degree of freedom in the calcu- 



lations (O'Regan et al. 2010), provided their localized, 



atomic character is maintained. It should also be pointed 
out that the value of U entering the calculations should 
not be considered universal, as it depends strongly on 
the manifold chosen (e.g. for pseudopotential calcula- 
tions on the oxidation state of the reference atomic cal- 



culation (Kulik and Marzari 2008)), or on the structure 
or electronic-structure of the problem studied. 

Last, it should be pointed out that functionals that 
aim to correct directly for some effects of self-interaction 



such as the Perdew-Zunger correction (Perdew and 



Zunger 1981 1 or Koopmans-compliant functionals (Dabo 



et al. 2010 ) - can lead naturally in a periodic system to 



Wannier-like localized orbitals that minimize the total 



energy (Park et al. 2011 Stengel and Spaldin 2008), 



while the canonical orbitals that diagonalize the Hamil- 
tonian still preserve Bloch periodicity. 



■51 



VIII. WANNIER FUNCTIONS IN OTHER CONTEXTS 



A. Phonons 



As described in Sec. |II.A{ Wannier functions provide an 
alternative, localized, description of a manifold of states 
spanned by the eigenstates (energy bands) of a single- 
particle Hamiltonian that describes electrons in a peri- 
odic potential. The equivalence of the Wannier repre- 
sentation and the eigenstate representation may be ex- 
pressed in terms of the band projection operator P, see 
Eq. (15). This operator satisfies the idempotency con- 
dition P 2 = P, which embodies simultaneously the re- 
quirements of orthogonality and Pauli exclusion. 

From their conception, and until relatively recently, 
Wannier functions have been used almost exclusively in 
this context, namely to represent a manifold of single- 
particle orbitals for electrons. Furthermore, as discussed 
in Sec. |nj we need not restrict ourselves to an isolated 
group of states, such as the occupied manifold: the dis- 
entanglement procedure enables a subspace of a larger 
manifold, e.g., of occupied and unoccupied states, to be 
selected which may then be wannierized. This has, for 
example, opened up areas of application in which Wan- 
nier functions are used as tight-binding basis functions 
for electronic structure and transport calculations, as de- 
scribed in Sec. Efl and Sec. NH\ 

From a general mathematical point of view, however, 
the set of orthogonal eigenfunctions of any self-adjoint 
(Hcrmitian) operator may be "rotated" by unitary trans- 
formation to another orthogonal basis that spans the 
same space. As we have seen, the unitary transforma- 
tion is arbitrary and may be chosen to render the new 
basis set maximally-localized, which has computational 
advantages when it comes to representing physical quan- 
tities that are short-ranged. When the operator in ques- 
tion has translational symmetry, the maximally-localized 
functions thus obtained are reminiscent of the Wannier 
functions familiar from electronic structure theory. Of- 
ten, such a basis is also preferable to using another local- 
ized basis because information regarding the symmetries 
of the self-adjoint operator from which the basis is de- 
rived is encoded within it. 

These ideas have led to the appropriation of the MLWF 
formalism described in Sec. |TT] for contexts other than 
the description of electrons: the single-particle electronic 
Hamiltonian is replaced by another suitable periodic self- 
adjoint operator, and the Bloch eigenstates are replaced 
by the eigenfunctions of that operator, which may then 
be transformed to give a MLWF-like representation that 
may be used as an optimal and compact basis for the de- 
sired calculation, for example, to analyze the eigenmodes 
of the static dielectric matrix in ice and liquid water ( Lu 



et al. 2008 1 



Below we review the three most prominent of these 
applications, namely to the study of phonons, photonic 
crystals, and cold atom lattices. 



Lattice vibrations in periodic crystals are usually de- 
scribed in terms of normal modes, which constitute a 
delocalized orthonormal basis for the space of vibra- 
tions of the lattice such that an arbitrary displace- 
ment of the atoms in the crystal may be expressed 
in terms of a linear combination of its normal modes. 



By analogy with the electronic case, Kohn (1973) first 



showed (for isolated phonon branches in one dimen- 
sion) that a similar approach could be used for con- 
structing a localized orthonormal basis for lattice vibra- 
tions that span the same space as the delocalized nor- 
mal modes. The approach was subsequently generalized 



to isolated manifolds in three-dimensions by Tindemans- 



van Eijndhoven and Kroese (1975). The localized modes 



are now gener ally referred to as lattice Wannier func 



tions (LWFs) (Iniguez et al 



1995) 



2000 



Rabe and Waghmare 



Following the notation of Sec. VI. D we denote by q the 
phonon wavevector, and by e q the matrix whose columns 
are the eigenvectors of the dynamical matrix. As in case 
of electronic Wannier functions, the phases of these eigen- 
vectors are undetermined. A unitary transformation of 
the form 



f.tV 



(134) 



performed within a subspace of dispersion branches that 
is invariant with respect to the space group of the crys- 
tal, results in an equivalent representation of generalized 
extended modes led that are also orthonormal. LWFs 
may then be defined by 



1 



iq R r~ 



with the associated inverse transform 



]T e iq R 

R 



1 1.1/ 



(135) 



(136) 



By construction, the LWFs are periodic according to 
% + t = wr, where t is a translation vector of the Born- 
von Karman supercell. 



The freedom inherent in Eq. ( 134 1 allows very local- 



ized LWFs to be constructed, by suitable choice of the 



transformation matrix Af q . As noted by Kohn (19731, 



the proof of exponential localization of LWFs follows 
the same reasoning as for electronic Wannier functions 



(SecpLG). 

The formal existence of LWFs was first invoked in order 
to justify the construction of approximate so-called local 
modes of vibration which were used in effective Hamilto- 
nians for the study of systems exhibiting strong coupling 
between electronic states and lattice instabilities, such as 



and Muller 1968) 



perovskite ferroelectrics (Pytte and Feder 1969 Thomas 
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Zhong et al. ( 1994 1 used first-principles methods in or- 



der to calculate the eigenvector associated with a soft 
mode at q = in BaTi03. A localized displacement 
pattern, or local mode, of the atoms in the cell was then 
parametrized, taking account of the symmetries associ- 
ated with the soft mode, and the parameters were fit- 
ted to reproduce the calculated soft mode eigenvector at 
q = 0. The degree of localization of the local mode was 
determined by setting to zero all displacement parame- 
ters beyond the second shell of atoms surrounding the 
central atom. Although this spatial truncation results in 
the local modes being non-orthogonal, it does not hamper 
the accuracy of practical calculations. As the local modes 
are constructed using information only from the eigen- 
vector at q = 0, they do not correspond to a particular 
phonon branch in the Brillouin zone. |Rabe and Wagh^| 
( 1995 ) generalized the approach to allow fitting to 



mare 



more than just q = 0, but rather to a small set of, usu- 
ally high-symmetry q-points. The phase-indeterminacy 
of the eigenvectors is exploited in order to achieve op- 
timally rapid decay o f the local modes. A nother ap- 



proach, introduced by Iniguez et al. (2000), constructs 



local modes via a projection method that preserves the 
correct symmetry. The procedure is initiated from sim- 
ple atomic displacements as trial functions. The quality 
of the local modes thus obtained may be improved by 
systematically densifying the q-point mesh that is used 
in Eq. (135). Although there is no formal criterion of 



maximal-localization in the approach, it also results in 
non-orthogonal local modes that decay exponentially. 

These ideas for generating local modes from first- 
principles calculations have been particularly success- 
ful for the study of structural phase transitions in 
ferroelectrics such as BaTiC>3 (Zhong et al 
PbTi0 3 (Waghmare and Rabe 



1995), 



Ni03 (Krakauer et al. 1999), Pb 3 GeTe 4 (Cockayne and 



Rabe 



2008) 



1994 



1997), Kb- 



1997) and perovskite superlattices (Lee et al. 



The use of maximal-localization as an exclusive cri- 
terion for determining LWFs was first introduced by 
Giustino and Pasquarello (2006). In this work, a 



real-space periodic position operator for non-interacting 
phonons was defined, by analogy with the periodic posi- 



tion operator for non-interacting electrons (Eq. (43)). 

The problem of minimizing the total spread of a set of 
WFs in real-space is equivalent to the problem of simulta- 
neously diagonalizing the three non-commuting matrices 
corresponding to the three components of the position 
operator represented in the WF basis, and |Giustino and 
Pasquarello ( 2006 1 use the method outlined by Gygi et al. 



(2003) to achieve this. It is worth noting that Giustino 



and Pasquarello ( 2006 1 furthermore define a generalized 



spread functional that, with a single parameter, allows 
a trade-off between localization in energy (the eigenstate 
or Bloch limit) and localization in space (the Wannier 
limit), resulting in so-called mixed Wannier-Bloch func- 



tions which may be obtained for the electrons as well as 
the phonons. 



Finally, as first pointed out by Kohn (1973), and sub- 



sequently by Giustino et al. (2007a), maximally-localized 



lattice Wannier functions correspond to displacements of 
individual atoms. This may be seen by considering a vi- 
brational eigenmode, e qs = e q;5 e lq R , and noting that it 
may be expressed as 



= J2 eiqR ' W<W< 



(137) 



s'R' 



Eq. (137) stands in direct correspondence to the elec- 



tronic analogue given by the inverse of Eq. (10), from 



which we conclude that the LWFs do indeed correspond 
to individual atomic displacements i5rr'(5 SS ' and, further- 
more, that the required unitary transformation is the ma- 
trix of eigenvectors [eq]^. As discussed in Sec. 



VI.D 



Giustino et al. (2007a) exploit this property for trie effi 



cient interpolation of dynamical matrices and calculation 
of electron-phonon couplings. 



B. Photonic crystals 

Photonic crystals are periodic arrangements of dielec- 
tric materials that are designed and fabricated in order to 



control the flow of light (John 1987 Yablonovitch 19871. 



They are very much to light what semiconductors are to 
electrons and, like semiconductors that exhibit an elec- 
tronic band gap in which an electron may not propagate, 
photonic crystals can be engineered to exhibit photonic 
band gaps: ranges of frequencies in which light is forbid- 
den to propagate in the crystal. In the electronic case, 
a band gap results from scattering from the periodic po- 
tential due to the ions in the crystal; in the photonic 
case, it arises from scattering from the periodic dielectric 
interfaces of the crystal. Again by analogy with elec- 
tronic materials, localized defect states can arise in the 
gap by the deliberate introduction of defects into a per- 
fect photonic crystal structure. The ability to control the 
nature of these states promises to lead to entirely light- 
based integrated circuits, which would have a number of 
advantages over their electronic counterparts, including 
greater speeds of propagation, greater bandwidth, and 



smaller energy losses ( Joannopoulos et al. 1997) 



In SI units, Maxwell's equations in source-free regions 
of space are 



V-E = 0, 



V ■ B = 0, 



<9B <9D 

VxE = - — , VxH = 



(138) 
(139) 



dt ' dt ' 

where the constitutive relations between the fields are 



D = £ r £oE, B = flj-flgH. 



(140) 



Considering non-magnetic materials (fi T = 1) with a po- 
sition dependent dielectric constant e r (r), and fields that 
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vary with a sinusoidal dependence e~ lult , it is straightfor- 
ward to derive electromagnetic wave equations in terms 
of either the electric field E or the magnetic field H, 

V x (V x E(r)) = ^e r (r)E(r), (141) 
Vx(e r - 1 (r)VxH(r)) = ^H(r), (142) 

where c = (^oeo) -1 ^ 2 is the speed of light. 

For a perfect periodic dielectric structure, e r (r) = 
e r (r + R), where R is a lattice vector. Application of 
Bloch's theorem leads to solutions that are indexed by 
wavevector k, which may be chosen to lie in the first 
Brillouin zone, and a band index n. For example 

H„ k (r)=e* k > ik (r), 

where u nk (r) = u nlc (r + R) is the periodic part of 
the magnetic field Bloch function. The electromagnetic 
wave equations can be solved, and hence the Bloch func- 
tions obtained, by a number of methods including finite- 



mnon 



difference time domain (Tafiove and Hagness 2005 Yee 



1966 ), transfer matrix (Pendry 1996 Pendry and Mack 



1992), empirical tight-binding methods (Lidorikis 



et al. 1998 Yariv et al. 19991, and Galerkin techniques 



in which the field is expanded in a set of orthogonal ba- 
sis functions ( Mogilevtsev et al. 1999 ) . Within the latter 
class, use of a plane-wave basis set is particularly com- 



Ho et al. 


1990 


Johnson and Joannopoulos 


2001) 



and, therefore, the fields satisfy orthogonality relations 
given b}0 

J rfrH* k (r)-H n ^(r) = 5 nn ^(k-k') ) (143) 
J dr Cr (r) E^(r) • EW(r) = 5 nn ,6(k - k'). (144) 



Leung (1993) first suggested that transforming to a 
basis of Wannier functions localized in real space would 
be advantageous for computational efficiency, especially 
when dealing with defects in photonic crystals which, us- 
ing conventional methods, require very large supercells 
for convergence. Although of great formal importance for 
justifying the existence of a suitable localized basis, and 
hence the tight-binding approach, the non-uniqueness of 
the transformation between Bloch states and Wannier 
functions caused difficulties. As a result, early work was 



limited to the case of single, isolated bands (Konotop 



1997 Leung 1993) or composite bands in which the ma- 
trix elements Umn were treated as parameters to fit the 
tight-binding band structure to the plane- wave result. 



21 The notation A-B = ~Yli—\ AiBi, and denotes the scalar product 
of the vectors A and B, with Cartesian components {Ai} and 
{Bi}, respectively. 



The formalism for obtaining maximally-localized Wan- 
nier functions, however, removed this obstacle and sev- 
eral applications of MLWFs to calculating the proper- 
ties of photonic crystals have been reported since, in 



both two-dimensional ( Garcia-Martin et al. 2003 Jiao 



et al. 2006 Whittaker and Croucher 2003) and three- 



dimensional (Takeda et al. 2006) photonic crystal struc- 



tures, as well as for the case of entangled bands (Hermann 



et al. 2008 ) (see Busch et al. ( 2003 ) for an early review) . 



Typically one chooses to solve for either the electric 
field E nk or the magnetic field H nk . Once the Bloch 
states for the periodic crystal are obtained, a basis of 
magnetic or electric field Wannier functions may be con- 
structed using the usual definition, e.g., for the magnetic 
field 

WSW— / dke-* k ' R 5>i k )H mk (r), (145) 
satisfying orthogonality relations 



rlr w' H '" 
Ul VV nR 



W 



(H) 
n'R' 



J RR' 



(146) 



where the unitary transformation Umn is chosen in the 
same way described in Sec. [H] such that the sum of the 
quadratic spreads of the Wannier functions is minimized, 
i.e., such that the Wannier functions are maximally lo- 
calized. 

Concentrating on the magnetic field, it may be ex- 
panded in the basis of Wannier functions with some ex- 
pansion coefficients c„r, 



H(r)=^c nR W^, ) (r), 



(147) 



uR 



which on substitution into Eq. ( 142 ) gives the tight- 
binding representation of the wave-equation for the mag- 
netic field in the Wannier function basis. 

The utility of the approach becomes evident when con- 
sidering the presence of a defect in the dielectric lattice 
such that e r (r) — > e r (r) + 8e(r). The magnetic field wave 
equations become 

V x ([e^r) + A-^r)] V x H(r)) = ^H(r), (148) 
where 

-5e(r) 



A(r) 



(149) 



e r (r)[e r (r) + Se(r)]- 
Using the Wannier functions from the defect-free calcula- 



tion as a basis in which to expand H(r), as in Eq. (147), 
the wave equations may be written in matrix form, 



iRR' 



(^m? + B nn'~ J ' 'it' - —'„].!- 



n'R' 



(150) 
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where 



C. Cold atoms in optical lattices 



aRR' _ 

nn' 



and 



V 



(2tt) 3 



BZ 



^.(R-RO^Wt 



u 



(k) 

mn' ' 



(151) 



B*$' = J drA(r) [Vx W„ R (r 



[Vx W B , R ,(r)], 
(152) 

Due to the localization and compactness of the basis, 
these matrix equations may be solved efficiently to find 
frequencies of localized cavity modes, dispersion rela- 
tions for waveguides, and the transmission and reflec- 



tion properties of complex waveguide structures. Fig. 35 



for example, shows the photonic band structure for a 
three-dimensional photonic crystal structure with a two- 
dimensional defect. 



x/a=0 

□ □ □ □ □ 

n n n n n 



k: 




FIG. 35 (Color online) Photonic band structure (bottom) of 
the 3D Si woodpile structure intercalated with a 2D layer 
consisting of a square lattice of square rods (top left). Solid 
lines indicate the photonic band structure calculated by the 
plane-wave expansion (PWE) method, and black points indi- 
cate that reproduced by the MLWFs. Shaded regions indicate 
the photonic band structure of the woodpile projected onto 
the 2D k|| space. The square rods in the 2D layer are cho- 
sen to structurally match the woodpile. The thickness of the 
layer is 0.8a, where a is the lattice parameter of the woodpile 
structure. Top right: absolute value of the 17th MLWF of the 
magnetic field in the yz plane. Adapted from |Takeda et al.\ 
(2006). 



A good 70 years after Albert Einstein predicted that a 
system of non-interacting bosons would undergo a phase 
transition to a new state of matter in which there is 
macroscopic occupation of the lowest energy quantum 
state, major achievements in laser cooling and evapora- 
tive cooling of atoms enabled the first experimental real- 



izations of Bosc-Einstcin condensation ( Anderson et al. 



1995 Bradley et al. 1995 Davis et al. 1995) and the 



award of Nobel prizes in 1997 and 2001. Since then, the 
study of cold atoms trapped in optical lattices has flour- 
ished. For a reviews see Bloch et al. (2008) and Morsch 



and Oberthaler (2006) 



Ultracold atoms trapped in optical lattices provide a 
versatile alternative to electrons in crystal lattices for the 
study of quantum phenomena. Indeed, they have a num- 
ber of advantages over the solid state in this respect, such 
as the absence of lattice defects, the absence of a counter- 
part to the electron-phonon interaction, and the ability 
to control precisely both the nature of the inter-atomic 
interactions and the depth and periodicity of the optical 
lattice potential. 

The second quantized Hamiltonian for a system of N 
weakly interacting bosons of zero spin and mass m in a 
(periodic) external potential Vb(r) = Vb(r + R) is given 
by ( |Yukalov[ [2009] ) 



U 



dv * t (i 



2m 



V 2 + F (r) 



*(r) 



dr * t (r)* t (r)#(r)#(r), 



(153) 



where g — 4ira s h 2 /m, it is assumed that the atoms 
interact via a short-range pseudopotential with a s as 
the s-wave scattering length, and ^(r) and 4d(r) are 
bosonic field operators obeying canonical commutation 



relations (Fetter and Walecka 2003). 



In a Bose-Einstein condensate, wherein the condensate 
particle densities are typically of the order of 10 14 cm~ 3 
or more, the mean-field limit of this Hamiltonian is usu- 
ally taken, which leads to the Gross-Pitaevskii equa- 
tion, also known as the non-linear Schrodinger equation 
(NLSE), 



^1 
2m 



V (r)+g\<p(T,t)\' 



<p{r,t). 



(154) 

yj(r, t) is the condensate wavefunction, the squared norm 
of which gives the condensate density. The Gross- 
Pitaevskii equation has been used with remarkable suc- 



cess in the study of BEC (Leggett 2001). 



As shown, for example, by Alfimov et al. (2002), a 



basis of Wannier functions, localized to each site a of the 
optical lattice, may be used to expand the condensate 
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wavef unction, 



(155) 



n,a 



The Wannier functions are related to the Bloch eigen- 
states ip n k{x) of the eigenvalue equation 



dx* +V ^> 



i>nk(x) = tnk^nk{x) (156) 



by the usual Wannier transformation 

w na (x) = ^[ dke l ^e- lkaL 4> nk {x). (157) 
^ Jbz 



Substituting Eq. (155) into Eq. (154), leads to a tight 



binding formulation known as the discrete non-linear 
Schrodinger equation (DNLSE), 



d_ 
dt 

where 



2tt 



. — i Al. p~ikuL 



(158) 
(159) 



BZ 



and the interaction matrix is given by 



UnijT = 9 I dxw na (x)w l p(x)w J1 (x)w kri (x). (160) 



Truncating the first term on the right-hand side of 



Eq. (158) to nearest-neighbors only, and the second term 
to on-site [a = ft = 7 = rj) terms within a single band 
(n = i = j = k) results in the usual tight-binding de- 



2001), 



scription ( Chiofalo et al. 20001 Trombettoni and Smerzi 



1+Cn,a+1 



(161) 

As pointed out by Alfimov et al. (2002), using a WF 



basis enables the range and type of interactions encapsu- 
lated in the DNLSE to be systematically controlled and 
improved. On the most part, however, WFs have been 
used in the context of the NLSE in order to carry out 
formal derivations and to justify the use of empirical or 
semi-empirical tight-binding models. 

An interesting analogy with electrons in atomic lat- 
tices manifests itself when the filling of sites in the optical 
lattice is low and hence particle correlations need to be 
accounted for more rigorously. This is done via the Bose- 



Hubbard model, developed by Fisher et al. ( 1989 ) in the 



context of He-4, and first applied to cold atoms in opti- 
cal lattices by Jaksch et al. (19981. The Bose- Hubbard 



Hamiltonian is derived from Eq. ( 153 ) by expanding the 



boson field operator in terms of WFs of a single band, 
localized at the lattice sites, 



(162) 



where the bosonic particle creation and annihilation op- 
erators, b a and W a , respectively, satisfy the usual com- 



mutation rules (Fetter and Walecka 2003). This, on 



approximation to nearest-neighbor coupling, and on-site 
only interactions, results in the standard Bose-Hubbard 



Hamiltonian (Jaksch et al. 1998) 



^bh = - J 2~Z Qfip + — ^a{n a - 1) 

<a,/3> a 



(163) 



where n a = b\J) a is the number operator for lattice site a, 
and the nearest-neighbor hopping and on-site repulsion 
parameters are given by 



J = — J dv wo(r) 



2m 



Vo(r) 



and 



U 



dr \w(r) 



wi(r), (164) 



(165) 



which may be calculated explicitly using WFs con- 



structed from Bloch eigenstates (Shotter et al. 2008 



Vaucher et al. 2007). The Bose-Hubbard model is the 
bosonic analogue to the Hubbard model for fermions. As 
in the latter case, the behavior of the model depends 
on the competition between hopping (J) and on-site (U) 
energies which determines whether the system is in a su- 
perfluid or a Mott insulator phase. 

Finally, we note that in work that is closely related 
to, and combines elements from both ideas developed in 
photonic crystals and cold atoms, WFs have also been 
used to represent polaritons in coupled cavity arrays, a 
class of systems that serves as another experimental re- 



alization of the Bose-Hubbard model (Hartmann et al. 



2008 2006) 



IX. SUMMARY AND PROSPECTS 

In this review, we have summarized methods for con- 
structing WFs to represent electrons in periodic solids 
or other extended systems. While several methods have 
been surveyed, our emphasis has been on the one of 



Marzari and Vanderbilt (1997), essentially the general 



ization of the approach of Foster and Boys (Boys 1960 



1966 Foster and Boys 1960a|b ) to periodic systems, 



in which the gauge freedom is resolved by minimizing 
the sum of the quadratic spreads of the WFs. The 
widespread adoption of this approach is reflected in the 
fact that it has been incorporated as a feature into a large 



5G 



number of modern first-principles electronic-structure 



code packages including Quantum ESPRESSO ( Gian- 



nozzi et al. 2009), Abinit (Gonze et al. 20091, FLEUR 



(Freimuth et al. 2008), WIEN2K (Kunes et al. 



Schwarz et al, 2002), Siesta (Korytar et al. 



Solcr et al. 2002), and VASP (Franchini et al 



Kresse and Furthmiiller 1996). In the above cases this 





2010 




2010 




2011 



has been done by providing an interface to the Wan- 



nier.90 package (Mostofi et al. 2008), an open-source 



post-processing code developed by the Authors, offering 
most of the capabilities described in this review. Other 
efforts have also seen the implementation of MLWFs in 



CPMD (CPMD 1990), GPAW (Enkovaara et al. 2010), 



OpenMX (|OpenMX[|2011[ ), and WanT ( |Calzolari et al. 
2004||Ferretti et aL||200"5af - this latter, and Wannier90, 
also allowing for quantum-transport calculations. 

After an initial wave of applications increased the vis- 
ibility of WFs in the community and demonstrated their 
utility for a variety of applications, other methods for 
constructing WFs were also developed, as discussed in 



Sees. |n| and III For some purposes, e.g., for many plane- 
wave based LDA+U and DMFT calculations, methods 
based on simple projection onto trial orbitals proved suf- 
ficient. Methods tuned specifically to T-point sampling 
of the BZ for supercell calculations also became popular. 
And, as surveyed briefly in Sec. |VIII| the construction 
and application of WFs was also extended to periodic 
systems outside the electronic-structure context, e.g., to 
phonons, photonic crystals, and cold-atom optical lat- 
tices. 

Still, the vast majority of applications of WF methods 
have been to electronic structure problems. The breadth 
of such applications can be appreciated by reviewing the 
topics covered in Sees. [TV|VII| Very broadly, these fall 
into three categories: investigations into the nature of 
chemical bonding (and, in complex systems such as liq- 
uids, the statistics of chemical bonding) , as discussed in 
Sec. |IV| applications that take advantage of the natural 
ability of WFs and WF charge centers to describe dipolar 
and orbital magnetization phenomena in dielectric, fer- 
roelectric, magnetic, and magnetoelectric materials, as 
reviewed in Sec. |Vj and the use of WF for basis func- 
tions, as surveyed in Sees. |VI|VII| 

Today these methods find applications in many topical 
areas including investigations into novel superconductors, 
multiferroics, and topological insulators. The importance 
of WFs is likely to grow in response to future trends in 
computing, which are clearly moving in the direction of 
more massive parallelization based on algorithms that 
can take advantage of real-space partitioning. This fea- 
ture of WFs should also facilitate their adoption in formu- 
lating new beyond-DFT methods in which many-body in- 
teractions are included in a real-space framework. Thus, 
the growing pressures for increased efficiency and accu- 
racy are likely to elevate the importance of WF-based 
methods in coming years. Overall, one can look forward 



to continued innovation in the development and applica- 
tion of WF-based methods to a wide variety of problems 
in modern condensed-matter theory. 
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